Activity in Pectoral Sandpipers: Statistics and figures

Background

Purpose: This notebook is used to perform statistical analyses and create tables/figures that are included in the manuscript.

Inputs: The following input datasets are required:

  • Capture data (data/capture-data.csv)
  • GPS data (data/gps-data.csv)
  • Nest data (data/nest-data.csv)
  • Parentage data (data/parentage-data.csv)
  • Sighting data (data/sighting-data.csv)
  • Cleaned ODBA data with activity thresholds (outputs/processed-data/ODBA-and-activity_10min-intervals_onboard.csv)
  • Male trips offsite (outputs/processed-data/male-trips-offsite.csv)

Data processing steps are outlined in the next section.

Outputs: Tables and figures for the manuscript are saved in “outputs/figures-and-tables”

Notes on data processing steps

Prior to performing the analyses below, steps for cleaning and preparing the data are performed in separate scripts:

  1. clean-onboard-odba-data.R: Cleans data processed on board the Druid accelerometer.
  2. define-activity-thresholds.Rmd: Determines activity thresholds for each individual
  3. map-male-locations.Rmd: Identifies when males took trips outside the study area

These processing scripts are found in the “scripts” folder. The key outputs are stored in “outputs/processed-data”.

Import data and prepare workspace

Most packages used in this Notebook can be installed using the install.packages function. One exception is “wadeR” which was developed by researchers at MPIO and needs to be installed from github.

# library(remotes)
# remotes::install_github("mpio-be/wadeR") # Utqiagvik (Barrow) mapping
library(tidyr) # for reshaping datasets
library(dplyr) # for other data manipulation
library(ggplot2) # for plotting
library(scales) # for plotting dates
library(lubridate) # for working with dates and times
library(ggalt) # for dumbbell charts
library(data.table) # for fast table reading and rolling functions
library(flextable) # generates neat tables
library(officer) # formatting flextable for Microsoft Word
library(fasttime) # for faster conversions to posixct
library(wadeR)
library(stringr) # count matches in string
library(here) # for easier directory navigation
library(glue) # for easier naming and referring to directories
library(dichromat) # for checking colours in figures
library(emmeans) # for computing estimated marginal means (EMMs)
library(effects) # for visualising effects
library(glmmTMB) # for glmm with beta error distribution
library(geosphere) # for distances lat lon
library(lme4) # for mixed effects models
library(spatialrisk) # for distances lat lon
library(patchwork) # better than ggarrange for arranging figures
library(DHARMa) # for checking model assumptions
library(sf) # for geometry functions
library(broom.mixed) # for converting glmmtmb to flextable

Some functions are used in multiple scripts and are therefore stored in their own script for easier access (general-helper-functions.R). This helper script needs to be loaded for the analyses below to run.

source("./scripts/general-helper-functions.R")

Male activity was recorded using miniature accelerometers and received in two different ways: processed on board the device (“onboard” data) or stored in raw form (“rawbased” data). During the course of the study, we also considered different methods for processing the data. These included summarising the data at either 1-min or 10-min intervals, and either defining activity thresholds for each individual separately (“individual” thresholds) or together (“global” thresholds).

There are advantages and disadvantages of different approaches, but results presented in the main text of the manuscript are based on the “onboard” data with a window-length of 10 minutes and “individual” thresholds.

In the supplementary material of the manuscript, we present a comparison of these different approaches. By changing the definitions in the code chunk below, and running the appropriate pre-processing steps, analyses in this R Notebook can be performed using datasets processed in different ways. These alternative datasets and pre-processing scripts are available upon request.

window_length = 10 # 1 or 10
data_type = "onboard" # rawbased or onboard
threshold_type = "individual" # global or individual

if(!(window_length==1 | window_length==10)){
  print("Error: Select a window length of either 1 or 10")
}

if(!(data_type=="onboard" | data_type=="rawbased")){
  print("Error: Select a data_type of either 'onboard' or 'rawbased'")
}

if(!(threshold_type=="individual" | threshold_type=="global")){
  print("Error: Select a threshold type of either 'individual' or 'global'")
}

if(threshold_type=="global" & window_length!=10){
  print("Warning: We recommend using a window length of 10 when using the global threshold. The global threshold was developed using data summarised at 10-minute intervals.")
}

if(threshold_type=="global" & data_type!="onboard"){
  print("Warning: We recommend using the 'onboard' data when applying the global threshold. The global threshold was developed using data processed on board the device.")
}

if(window_length=="1" & data_type!="onboard"){
  print("Error: Onboard data is only available with a window length of 10 minutes.")
}

if(threshold_type == "global") {
  global_threshold = 0.27 # this threshold is derived from a separate analysis that is available upon request
}
# 10-min onboard activity data is used for the main analyses in the manuscript
# all other data types are stored in a supplementary folder

if(window_length==10 & data_type=="onboard") {
  activity_data_path = here("outputs","processed-data")
} else {
  activity_data_path = here("supplementary", "supp-data")
}
captures <- read.csv(here("data","capture-data.csv"))
nests <- read.csv(here("data","nest-data.csv"))
sightings <- read.csv(here("data", "sighting-data.csv"))
parentage <- read.csv(here("data", "parentage-data.csv"))
gps <- read.csv(here("data","gps-data.csv"))

activity_onboard <- read.csv(here("outputs", "processed-data",
                                  "ODBA-and-activity_10min-intervals_onboard.csv"))
  # this dataset is used to help determine when final data transmission occurred

activity <- read.csv(glue("{activity_data_path}/ODBA-and-activity_{window_length}min-intervals_{data_type}.csv"))
  # this dataset is the source of ODBA and activity data for analyses
  # there is the option to customise the input to analyse different activity datasets

male_trips_offsite <- read.csv(here("outputs", "processed-data", "male-trips-offsite.csv"))

All dates and times are local Alaskan times, unless specified otherwise.

captures <- captures %>%
  mutate(date=local_date_to_posix(as.character(date)),
         released_dt = local_dt_to_posix(paste(date,released)))

nests <- nests %>%
  mutate(date_found=local_date_to_posix(as.character(date_found)),
         est_hatch_date=local_date_to_posix(as.character(est_hatch_date)),
         hatch_date=local_date_to_posix(as.character(hatch_date)))

sightings <- sightings %>%
  mutate(datetime=local_dt_to_posix(as.character(datetime)))

gps <- gps %>%
  mutate(datetime=local_dt_to_posix(as.character(datetime)))

activity_onboard <- activity_onboard %>%
  mutate(window_end_utc=acc_to_posix(as.character(window_end_utc)),
         window_end_localtime=dt_to_local(window_end_utc))

activity <- activity %>%
  mutate(window_end_utc=acc_to_posix(as.character(window_end_utc)),
         window_end_localtime=dt_to_local(window_end_utc))

male_trips_offsite <- male_trips_offsite %>%
  mutate(start_dt_local=local_dt_to_posix(as.character(start_dt_local)),
         end_dt_local=local_dt_to_posix(as.character(end_dt_local)))

To avoid biases that might arise from incomplete data, we use a minimum threshold for data availability for some analyses. Specifically, we have defined a minimum number of hours of data that need to be available per individual per day (“min_acc_hrs_per_day”) and a minimum number of days of data that need to be available per individual (“min_acc_days”). Data that do not meet these requirements (i.e. insufficient data available from an individual on a particular day, or insufficient data available from a particular individual overall) are excluded from some analyses. By default, we have required at least 85% complete data per individual per day, and at least 2 days of data from each individual.

min_acc_hrs_per_day = 20.4 # 20.4 means at least 85% complete
min_acc_days = 2 # except where otherwise specified
male_colour = "darkgreen"
female_colour = "gray32" 
success_colour = "#2267E7"  #"#20EAABFF"
fail_colour = "#DF536B"
default_shape = 21
uncertain_shape = 4

legend.heading.size = 10
axis.heading.size = 10
label.size = 10
font.family = "Calibri"
small.point.size = 1

# ribbon plot specifications
point_size = small.point.size+1
point_alpha = 0.5
col_scale <- c("springgreen4", "darkgray")
extra_nest_shape = 13

windowsFonts("Calibri" = windowsFont("Calibri")) # avoids font issues
set_flextable_defaults(
  font.size = 11,
  theme_fun = theme_booktabs,
  font.family = "Calibri",
  padding = 3,
  background.color = "white")

sect_properties <- prop_section(
  page_size = page_size(orient = "portrait",
                        width = 21 * 0.393701,
                        height =29.7  * 0.393701),
  type = "continuous",
  page_margins = page_mar(bottom = 2.54 * 0.393701,
                          top = 2.54 * 0.393701,
                          left = 2.54 * 0.393701,
                          right = 2.54 * 0.393701)
)
save_plots = TRUE
save_tables = TRUE
output_path <- here('outputs','figures-and-tables')
output_supp_effects <- here('supplementary','supp-outputs','effect-sizes')

if (!dir.exists(output_path)) {
  dir.create(output_path)
}

Preparing data for analysis

Note: The tables and figures in this section are for context and data-checking only. They are not included in the manuscript.

Identify tagged males

tagged_males <- captures %>%
  filter(tag_action=="on") %>% 
  filter(!colour_combo %in% colour_combos_to_exclude) %>%
  select(colour_combo, bird_ID, tag_ID, date, released_dt, body_mass)
ggplot(tagged_males, aes(x=date))+
  theme_classic()+
  geom_bar(stat="count",
           fill=male_colour) +
  scale_x_datetime(breaks=date_breaks(width="1 day"),
               date_labels="%d",
               expand = c(0, 0)) +
  coord_cartesian(xlim=c(local_date_to_posix("2022-06-06"),
                         local_date_to_posix("2022-06-28")))+
  xlab("\n Date (June)")+
  ylab("Number of new tags deployed on males\n")

In 2022, we equipped 100 adult male pectoral sandpipers with 3.7 g accelerometer tags (weight includes harness and glue) with GPS function (Debut NANO, Druid Inc.) between 2022-06-07 and 2022-06-20. We targeted males who were exhibiting territorial behaviours. We stopped tagging males when we ran out of tags (although new males were still entering the study site).

The weight of the device was between approximately 3.2 and 4.2 % of the bird’s body mass.

Determing female fertility

We defined the fertile period of each female pectoral sandpiper as the interval beginning 5 days before the female began laying eggs (Kempenaers and Valcu 2017 Nature, https://doi.org/10.1038/nature20813) and ending the day the penultimate egg was laid. Clutch size ranged from 2 to 4 eggs (89.7% of nests contained 4 eggs, 6.9% contained 3 eggs, and 3.4% contained 2 eggs), meaning that the fertile period for each female ranged from 6-8 days. It is not actually known whether females still engage in copulations after initiating their clutch.

If the nest was found during laying (clutch size when found < maximum clutch size), the beginning of the fertile period was defined as date found - clutch size when found (females lay one egg per day).

If the nest was complete when found, the beginning of the fertile period was defined as: hatch date - the typical incubation period (22 days) - (maximum clutch size - 1) . For hatch date, we used either the earliest actual hatch date if the eggs successfully hatched in the incubator, or the earliest predicted hatch date (based on the floatation method) if they did not.

Our calculations assume that whenever the nest was found or checked, the female had not laid any eggs yet on that day.

nests$found_status <- if_else(
  nests$first_clutchsize==nests$max_clutchsize, 
  "complete",
  "laying"
)

nests <- nests %>%
  mutate(clutch_init = case_when(found_status=="laying" ~ date_found-first_clutchsize*(24*60*60),
                                 is.na(hatch_date) ~ est_hatch_date - 22*(24*60*60) - (max_clutchsize-1)*(24*60*60),
                                 TRUE ~ hatch_date - 22*(24*60*60) - (max_clutchsize-1)*(24*60*60)),
         penult_egg = case_when(max_clutchsize>1 ~ clutch_init + (max_clutchsize-2)*(24*60*60),
                                TRUE ~ clutch_init),
         first_fertile_day=clutch_init-5*(24*60*60),
         last_fertile_day = penult_egg)

For the overall study site, we defined the start of the fertile period as the 5% quantile of all fertile start dates, and the end as the 95% quantile. A similar approach was used to define the end of the fertile period in Lesku et al. 2012 Science (https://doi.org/10.1126/science.1220939.

After the start date is ‘fertile’, end date onwards is ‘post-fertile’.

start_FertilePeriod <- floor_date(
  quantile(nests$first_fertile_day, 0.05, na.rm=T),
  unit="day")
end_FertilePeriod <- floor_date(
  quantile(nests$last_fertile_day, 0.95, na.rm=T),
  unit="day")
all_FertileDates <- seq(min(nests$first_fertile_day, na.rm=T), 
                        by = "day", 
                        length.out = max(nests$last_fertile_day, na.rm=T)-
                          min(nests$first_fertile_day, na.rm=T))

countFertileFemales <- function(date) {
  count <- date>=nests$first_fertile_day &
    date<=nests$last_fertile_day
  sum <- sum(count, na.rm=TRUE)
  return(sum)
}

ff_by_date <- data.frame(date=unlist(all_FertileDates))
ff_by_date$n_FertileFemales <- sapply(all_FertileDates,
                                      countFertileFemales)

barplot_ff <- ggplot(ff_by_date, aes(x=date, 
                                     y=n_FertileFemales)) +
  # theme_classic()+
  theme(panel.background=element_blank(),
        axis.line = element_line(size=0.5, colour="black"))+
  annotate(geom="rect",
           xmin=start_FertilePeriod-24*60*30,
           xmax=end_FertilePeriod+24*60*30,
           ymin=0,
           ymax=43,
           fill="grey",
           alpha=0.3)+
  geom_bar(stat="identity",
           fill=female_colour) +
  coord_cartesian(xlim=c(local_date_to_posix("2022-06-06"),
                         local_date_to_posix("2022-07-02")))+
  scale_y_continuous(limits=c(0,43), expand=c(0,0))+
  xlab("\n Date") +
  ylab("Number of fertile females\n")

barplot_ff

Summarise paternity data

We collected blood samples from all of the males we tagged (as well as some that we didn’t) and from eggs/chicks on our study plot. DNA was extracted from these samples by collaborators at the University of Alaska. These DNA samples were shipped to the research institute.

parentage <- left_join(parentage, captures[,c("bird_ID","colour_combo")], 
                       by=join_by("ID.father"=="bird_ID")) %>%
  mutate(ID.father=as.character(ID.father)) %>%
  mutate(ID.father=replace_na(ID.father,"unbanded")) %>%
  rename(LL.father=colour_combo) %>%
  mutate(LL.father=replace_na(LL.father,"unbanded")) %>%
  left_join(., captures[,c("bird_ID","colour_combo")], by=join_by("ID.mother"=="bird_ID")) %>%
  rename(LL.mother=colour_combo)
fathers <- parentage %>%
  group_by(ID.father, LL.father, Nest) %>%
  summarise(n_eggs=n()) %>%
  group_by(ID.father, LL.father) %>%
  summarise(n_eggs=sum(n_eggs),
            n_nests=n())

fathers %>% flextable()

ID.father

LL.father

n_eggs

n_nests

141269009

PI,PI,R

6

2

141269026

O,DB,O

4

1

141269043

DG,O,O

4

1

141269044

Y,DG,PI

5

2

141269046

DB,DG,Y

4

1

141269114

PI,R,PI

3

1

141269125

Y,Y,R

1

1

141269127

DG,DB,DG

4

1

141269144

DB,Y,R

3

1

141269146

R,DG,PI

2

1

141269147

O,DG,DG

5

2

141269148

DB,R,O

3

1

141269152

Y,Y,Y

1

1

141269169

DB,R,R

4

1

141269175

DG,R,R

3

1

141269182

O,Y,O

2

1

unbanded

unbanded

136

38

(Table for background only; not to be included in paper)

parentage %>% 
  group_by(Nest) %>% 
  summarise(fathers=toString(unique(LL.father)),
            unique_fathers=n_distinct(LL.father)) %>% 
  filter(fathers!="NA" & fathers!="unbanded") %>%
  arrange(desc(unique_fathers)) %>%
  flextable()

Nest

fathers

unique_fathers

P1103

unbanded, DG,R,R

2

P1401

unbanded, Y,Y,Y

2

P201

unbanded, Y,DG,PI

2

P208

unbanded, O,Y,O

2

P403

PI,R,PI, O,DG,DG

2

P701

R,DG,PI, PI,PI,R

2

P702

Y,Y,R, DB,Y,R

2

P1001

DG,DB,DG

1

P1003

O,DB,O

1

P1403

Y,DG,PI

1

P1903

DB,DG,Y

1

P1905

DG,O,O

1

P203

DB,R,O

1

P204

DB,R,R

1

P214

O,DG,DG

1

P401

PI,PI,R

1

Exclude unreliable or irrelevant GPS data

GPS data are filtered down to data collected while the device was on the bird.

Occasionally, the devices recorded a latitude and longitude of 200. This is an error value that doesn’t mean anything and should be excluded.

male_gps <- gps %>%
  left_join(., tagged_males, join_by('tag_ID')) %>%
  filter(lat!=200 & lon!=200) %>% # these are nonsense values
  filter(datetime>released_dt) %>%
  mutate(local_date=floor_date(datetime, unit="days"),
         bird_ID=as.factor(bird_ID)) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = FALSE) 

Determining site tenure

We use capture date as a proxy for arrival date.

We use the latest of sighting, GPS and ODBA data as a proxy for departure.

There are two birds who were “seen” after they last transmitted data. One of these was definitely a mistake; the wrong GPS point must have been recorded, because the sighting is associated with a video file from 10 days earlier (14 June rather than 24 June). So for this bird (DG,DB,DG) I am using the last date of data transmission as his departure date.

The other I have no reason to exclude, and the sighting was only 1 day after he last transmitted data.

latest_onboard = activity_onboard %>%
  group_by(tag_ID) %>%
  mutate(window_end_localtime_hour=floor_date(window_end_localtime, unit="hours")) %>%
  summarise(last_ODBA = max(window_end_localtime, na.rm=T),
            last_hour_n_ODBA = sum(window_end_localtime>=(last_ODBA-60*60), na.rm=T))

latest_gps = male_gps %>%
  group_by(tag_ID) %>%
  summarise(last_GPS = max(datetime, na.rm=T),
            last_hour_n_GPS = sum(datetime>=(last_GPS-60*60), na.rm=T))

latest_sightings <- inner_join(sightings, tagged_males, by='colour_combo') %>%
  group_by(tag_ID) %>%
  summarise(last_sighting = max(datetime, na.rm=T),
            last_hour_n_sightings = sum(datetime>=(last_sighting-60*60), na.rm=T))

male_departures = left_join(tagged_males, latest_onboard, join_by("tag_ID")) %>%
  left_join(., latest_gps, join_by("tag_ID")) %>%
  left_join(., latest_sightings, join_by("tag_ID"))

male_departures = male_departures %>%
  mutate(last_data_akdt = pmax(last_ODBA,last_GPS,last_sighting, na.rm=T),
         last_data_utc = dt_to_acc(last_data_akdt),
         last_data_cet = with_tz(last_data_akdt, tz="CET"),
         across(c(last_ODBA, last_GPS, last_sighting), 
                ~case_when(. == last_data_akdt ~ cur_column()), .names = 'new_{col}')) %>%
  unite(col=last_data_type, starts_with('new'), na.rm = TRUE, sep = ', ') %>%
  mutate(last_hour_n = case_when(last_data_type=="last_ODBA" ~ last_hour_n_ODBA,
                                    last_data_type=="last_GPS" ~ last_hour_n_GPS,
                                    last_data_type=="last_sighting" ~ last_hour_n_sightings,
                                    TRUE ~ NA))

male_tenure = male_departures %>%
  select(colour_combo, bird_ID, tag_ID, released_dt, last_data_akdt, last_data_type, last_hour_n) %>%
  rename(arrival_dt=released_dt,
         departure_dt=last_data_akdt) %>%
  mutate(arrival_date=floor_date(arrival_dt, unit="days"),
         departure_date=floor_date(departure_dt, unit="days"),
         departure_dt_confidence = case_when(last_hour_n>1 & last_data_type!="last_sighting" ~ "OK",
                                          TRUE ~ "low"),
         unrevised_tenure_days_exact = as.numeric(difftime(departure_dt, arrival_dt, unit="days")),
         unrevised_tenure_days = as.numeric(difftime(departure_date, arrival_date, unit="days")))

rm(male_departures, latest_onboard, latest_gps, latest_sightings)

Some males left and then returned again. We subtract the days that they were away from tenure.

male_trips_offsite <- male_trips_offsite %>%
  mutate(offsite_days = as.numeric(difftime(as.Date(end_dt_local),
                                  as.Date(start_dt_local),
                                  unit="days")))

male_trips_offsite %>%
  summarise(mean_offsite_days=mean(offsite_days),
            sd_offsite_days=sd(offsite_days),
            min_offsite_days=min(offsite_days),
            max_offsite_days=max(offsite_days))
##   mean_offsite_days sd_offsite_days min_offsite_days max_offsite_days
## 1              3.57            2.07                1                7
male_offsite_summary = male_trips_offsite %>%
  select(bird_ID, duration_hrs, offsite_days) %>%
  rename(offsite_hrs=duration_hrs)
male_tenure <- left_join(male_tenure, male_offsite_summary, join_by('bird_ID')) %>%
  mutate(tenure_days_exact=case_when(is.na(offsite_hrs) ~ unrevised_tenure_days_exact,
                                     TRUE ~ unrevised_tenure_days_exact-(offsite_hrs/24)),
         tenure_days=case_when(is.na(offsite_days)~unrevised_tenure_days,
                               TRUE ~ unrevised_tenure_days-offsite_days)) %>%
  select(colour_combo, bird_ID, tag_ID, 
         arrival_dt, departure_dt, arrival_date, departure_date,
         tenure_days, tenure_days_exact)

Examine breeding density relative to other years

In total, there were 58 breeding females and median male tenure was 10 days. This was therefore a moderate density breeding season, compared to other years (Kempenaers and Valcu 2017, Nature; number of breeding females ranges from approx. 23-90). It was higher density year than 2008 and 2009, which were the breeding seasons when data were collected for the previous activity study (Lesku et al. 2012, Science: 39 nests in 2008 and 17 nests in 2009.)

Exclude unreliable or irrelevant ODBA data

A simple and common measure of activity is overall dynamic body acceleration (ODBA). ODBA is a measure of total movement, combined across the three axes of the accelerometer. It is calculated by first subtracting static acceleration from each axis (acceleration related to posture; most often calculated by taking the 2-second running mean) to get dynamic acceleration, then summing together the absolute dynamic acceleration for each axis (for example, see: https://besjournals.onlinelibrary.wiley.com/doi/10.1111/j.2041-210X.2010.00057.x) ODBA is also commonly used as a proxy for energy expenditure related to movement.

Data have already been excluded when the tag was not on the bird. In a separate script ( define-activity-thresholds.RmD ) I also defined an activity threshold for each individual bird.

Below, I am excluding the following data:

  • Activity thresholds from birds without reliable estimates (due to <36 hours total available data)
  • Data collected when birds were offsite
  • Data collected <12 hours after putting the tag on
min_hrs_after_release = 12
# Exclude activity thresholds with <36 hours of data (when using individual threshold)
if(threshold_type=="individual") {
activity <- activity %>%
  mutate(activity_threshold=case_when(n_windows_threshold<=36*60/window_length ~ NA,
                                      TRUE ~ activity_threshold),
         active=case_when(mean_ODBA>activity_threshold ~ "active",
                          mean_ODBA<=activity_threshold ~ "inactive",
                          TRUE ~ NA))
}

# Exclude trips offsite
activity <- male_trips_offsite %>%
  select(colour_combo, bird_ID, start_dt_local, end_dt_local) %>%
  right_join(., activity, join_by('colour_combo','bird_ID'),
             relationship="one-to-many") %>%
  filter(is.na(start_dt_local) |
         window_end_localtime < start_dt_local |
         window_end_localtime-window_length*60 > end_dt_local) %>%
  select(!start_dt_local & !end_dt_local)

# Exclude 12 hours after release
activity <- tagged_males %>%
  select(colour_combo, bird_ID, released_dt) %>%
  right_join(., activity, join_by('colour_combo','bird_ID')) %>%
  filter(window_end_localtime-window_length*60 > released_dt+(min_hrs_after_release*60*60))

There was also one male who travelled 300 km in one day (out over the ocean for 150km, then back; see map-male-locations.Rmd). It doesn’t make sense to include this data in our analysis. He was R,Y,DB and his travel date was 20 June.

activity <- activity %>%
  mutate(window_end_local_date=floor_date(window_end_localtime, unit="day"))

activity <- activity %>%
  filter(!(colour_combo=="R,Y,DB" & window_end_local_date==as.Date("2022-06-20")))

Summarise and define activity thresholds and levels

if(threshold_type=="individual") {
activity %>% 
  summarise(min_threshold=min(activity_threshold, na.rm=T),
            max_threshold=max(activity_threshold, na.rm=T),
            median_threshold=median(activity_threshold, na.rm=T)) %>%
  flextable()
} else {
  print(paste0("Global threshold = ", global_threshold))
}

min_threshold

max_threshold

median_threshold

0.147

0.37

0.239

By determining an activity threshold for each individual, based on the distribution of his activity data, we should minimise bias in activity estimates caused by variation in measurements by different devices and minor variations in attachment. However, it is possible that we are also sometimes underestimating the activity of highly active birds and overestimating the activity of less active birds, since these distributions in activity distributions can affect the values of the threshold. Compared to having a global threshold across all birds, this method is therefore more conservative and more likely to underestimate inter-individual variation in activity levels.

Note: Birds without activity thresholds (i.e. <36 hours data) have not been excluded from the “activity” dataset; for some analyses, the ODBA values of these birds are still relevant and worth including.

if(threshold_type=="global") {
activity <- activity %>%
  mutate(activity_threshold=global_threshold,
         active=case_when(mean_ODBA>global_threshold ~ "active",
                          TRUE ~ "inactive"))
}

Calculate daily activity

activity$local_date = floor_date(activity$window_end_localtime, unit="days")

daily_activity <- activity %>%
  #filter(!is.na(activity_threshold)) %>%
  group_by(bird_ID, tag_ID, colour_combo, local_date) %>%
  summarise(activity_threshold=first(activity_threshold),
            n_windows=n(),
            n_hours=n_windows*window_length/60,
            perc_active=sum(active=="active")/n_windows*100,
            total_active=sum(active=="active"),
            mean_ODBA=mean(mean_ODBA))

Determine distance of males from fertile nests

predated_nests <- nests %>%
  filter(is.na(est_hatch_date)) %>% # predated nests had no estimated hatch date, eggs not floated
  select(nest, date_found, first_clutchsize, lat, lon)

predated_nests %>% flextable()

nest

date_found

first_clutchsize

lat

lon

P1109

2022-06-30 00:00:00

4

71.3

-157

P1110

2022-06-30 00:00:00

3

71.3

-157

P1111

2022-06-30 00:00:00

4

71.3

-157

P1112

2022-07-01 00:00:00

2

71.3

-157

P1409

2022-06-28 00:00:00

4

71.3

-157

P1910

2022-06-30 00:00:00

4

71.3

-157

P206

2022-06-30 00:00:00

4

71.3

-157

P207

2022-06-30 00:00:00

4

71.3

-157

To estimate the fertile dates of predated nests, I’ve assumed that the maximum clutch size would have been 4 (so nests found with clutch size < 4 were found during laying).

For nests that were not found during laying, we only know that the female was no longer fertile when the nest was found. All of the nests were found after the end of the fertile phase for the plot (based on the fertile periods of other nests). I am therefore using the entire fertile phase for the plot as the fertile period for these predated nests. This is used to help establish which males were spending time near the nest when it might have been fertile.

predated_nests <- predated_nests %>%
  mutate(clutch_init = case_when(first_clutchsize<4 ~ date_found-first_clutchsize*86400,
                                 TRUE ~ NA),
         penult_egg = case_when(is.na(clutch_init) ~ NA,
                                TRUE ~ clutch_init + 2*86400),
         first_fertile_day = case_when(!is.na(clutch_init) ~ clutch_init-5*86400,
                                       TRUE ~ start_FertilePeriod),
         last_fertile_day = case_when(!is.na(penult_egg) ~ penult_egg,
                                      TRUE ~ end_FertilePeriod)
  )

predated_nests %>% flextable()

nest

date_found

first_clutchsize

lat

lon

clutch_init

penult_egg

first_fertile_day

last_fertile_day

P1109

2022-06-30 00:00:00

4

71.3

-157

2022-06-10 00:00:00

2022-06-26 00:00:00

P1110

2022-06-30 00:00:00

3

71.3

-157

2022-06-27 00:00:00

2022-06-29 00:00:00

2022-06-22 00:00:00

2022-06-29 00:00:00

P1111

2022-06-30 00:00:00

4

71.3

-157

2022-06-10 00:00:00

2022-06-26 00:00:00

P1112

2022-07-01 00:00:00

2

71.3

-157

2022-06-29 00:00:00

2022-07-01 00:00:00

2022-06-24 00:00:00

2022-07-01 00:00:00

P1409

2022-06-28 00:00:00

4

71.3

-157

2022-06-10 00:00:00

2022-06-26 00:00:00

P1910

2022-06-30 00:00:00

4

71.3

-157

2022-06-10 00:00:00

2022-06-26 00:00:00

P206

2022-06-30 00:00:00

4

71.3

-157

2022-06-10 00:00:00

2022-06-26 00:00:00

P207

2022-06-30 00:00:00

4

71.3

-157

2022-06-10 00:00:00

2022-06-26 00:00:00

nests <- nests %>%
  filter(!nest %in% predated_nests$nest) %>%
  full_join(., predated_nests, 
            join_by('nest', 'date_found', 'lat','lon','first_clutchsize', 'clutch_init', 'penult_egg','first_fertile_day','last_fertile_day'))
nest_parentage <- parentage %>%
  filter(is.na(Comments) | Comments!="could not be genotyped") %>%
  select(Nest, LL.father) %>%
  unique() %>%
  group_by(Nest) %>%
  arrange(LL.father) %>%
  summarise(father1 = first(LL.father),
            father2 = case_when(last(LL.father)==father1 ~ NA,
                                TRUE ~ last(LL.father)),
            n_fathers = n()) %>%
  rename(nest = Nest)
nest_dates <- nests %>%
  select(nest, first_fertile_day, last_fertile_day) %>%
  pivot_longer(cols = first_fertile_day:last_fertile_day,
               values_to = 'date') %>%
  filter(!is.na(date)) %>%
  group_by(nest) %>%
  complete(date = seq(min(date), max(date), by='days')) %>%
  select(-name)

nest_long <- left_join(nest_dates, nests, join_by('nest'), relationship='many-to-one') %>%
  left_join(., nest_parentage, join_by('nest'), relationship='many-to-one') %>%
  select(nest, date, first_fertile_day, last_fertile_day, lat, lon, father1, father2) %>%
  st_as_sf(coords = c("lon", "lat"), crs = 4326, remove = FALSE)
male_date_gps <- male_gps %>%
  group_by(colour_combo, bird_ID, tag_ID, local_date, released_dt) %>%
  summarise(n_gps=n(),
            centroid = st_centroid(st_union(geometry)))
daily_activity_gps <- activity %>%
  group_by(bird_ID, tag_ID, colour_combo, local_date) %>%
  summarise(activity_threshold=first(activity_threshold),
            mean_ODBA=mean(mean_ODBA),
            perc_active=case_when(is.na(activity_threshold) ~ NA,
                                  TRUE ~ sum(active=="active")/(sum(active=="active")+sum(active=="inactive"))*100),
            prop_activity=case_when(is.na(activity_threshold) ~ NA,
                                    TRUE ~ perc_active/100),
            total_active_windows=sum(active=="active"),
            n_windows=n(),
            n_hours=n()/(60/window_length)) %>%
  mutate(bird_ID=as.factor(bird_ID)) %>%
  filter(n_hours>min_acc_hrs_per_day) %>%
  select(bird_ID, tag_ID, colour_combo, local_date, mean_ODBA, perc_active, prop_activity, n_windows) %>%
  full_join(., male_date_gps, join_by('tag_ID', 'bird_ID', 'colour_combo', 'local_date'), 
            relationship='many-to-many') %>%
  filter(!is.na(mean_ODBA))

Note that this next chunk of code takes a while to run (getting distance from each male from each nest on each day; ~ 3-5 mins).

nest_male_date <- inner_join(nest_long,
                             daily_activity_gps,
                             join_by('date'=='local_date'),
                             relationship='many-to-many') %>%
  filter(!is.na(lat) & !is.na(mean_ODBA))

nest_male_date <- nest_male_date %>%
  rowwise() %>%
  mutate(dist_m = st_distance(geometry, centroid))
nest_male_date <- nest_male_date %>%
  filter(!is.na(n_gps)) %>%
  mutate(success=case_when(colour_combo==father1 ~ 1,
                           colour_combo==father2 ~ 1,
                           TRUE ~ 0),
         bird_ID = as.factor(bird_ID),
         date = as.factor(date),
         nest = as.factor(nest))
nest_male_date %>%
  filter(success==1) %>%
  select(nest, father1, father2, dist_m, colour_combo, success) %>%
  group_by(nest, father1, father2, colour_combo, success) %>%
  summarise(mean_dist_m=mean(dist_m),
            min_dist_m=min(dist_m),
            median_dist_m=median(dist_m),
            max_dist_m=max(dist_m),
            n=n()) %>%
  arrange(median_dist_m) %>%
  select(!c(geometry)) %>%
  flextable()

nest

father1

father2

colour_combo

success

mean_dist_m

min_dist_m

median_dist_m

max_dist_m

n

geometry

P208

O,Y,O

unbanded

O,Y,O

1

34.6 [m]

8.7 [m]

33.3 [m]

52.7 [m]

7

[[XY]]

P203

DB,R,O

DB,R,O

1

38.8 [m]

19.9 [m]

36.0 [m]

73.6 [m]

7

[[XY]]

P214

O,DG,DG

O,DG,DG

1

273.9 [m]

31.5 [m]

57.7 [m]

1,035.8 [m]

8

[[XY]]

P1403

Y,DG,PI

Y,DG,PI

1

62.1 [m]

62.1 [m]

62.1 [m]

62.1 [m]

1

[[XY]]

P1903

DB,DG,Y

DB,DG,Y

1

84.2 [m]

47.8 [m]

63.0 [m]

163.1 [m]

4

[[XY]]

P1001

DG,DB,DG

DG,DB,DG

1

66.5 [m]

66.5 [m]

66.5 [m]

66.5 [m]

1

[[XY]]

P403

O,DG,DG

PI,R,PI

O,DG,DG

1

69.9 [m]

63.8 [m]

69.5 [m]

76.3 [m]

3

[[XY]]

P1003

O,DB,O

O,DB,O

1

91.4 [m]

79.4 [m]

92.1 [m]

102.8 [m]

7

[[XY]]

P204

DB,R,R

DB,R,R

1

112.5 [m]

105.5 [m]

111.1 [m]

120.8 [m]

3

[[XY]]

P201

Y,DG,PI

unbanded

Y,DG,PI

1

289.1 [m]

43.2 [m]

118.8 [m]

799.9 [m]

5

[[XY]]

P1905

DG,O,O

DG,O,O

1

165.3 [m]

112.1 [m]

119.1 [m]

383.9 [m]

6

[[XY]]

P401

PI,PI,R

PI,PI,R

1

172.5 [m]

142.0 [m]

162.8 [m]

214.6 [m]

7

[[XY]]

P702

DB,Y,R

Y,Y,R

Y,Y,R

1

173.7 [m]

148.6 [m]

173.7 [m]

198.7 [m]

2

[[XY]]

P1401

Y,Y,Y

unbanded

Y,Y,Y

1

348.1 [m]

167.0 [m]

187.8 [m]

1,350.2 [m]

8

[[XY]]

P403

O,DG,DG

PI,R,PI

PI,R,PI

1

191.0 [m]

135.2 [m]

208.0 [m]

222.9 [m]

7

[[XY]]

P701

PI,PI,R

R,DG,PI

PI,PI,R

1

763.5 [m]

741.1 [m]

768.0 [m]

777.1 [m]

4

[[XY]]

P1103

DG,R,R

unbanded

DG,R,R

1

800.7 [m]

800.7 [m]

800.7 [m]

800.7 [m]

1

[[XY]]

P701

PI,PI,R

R,DG,PI

R,DG,PI

1

896.3 [m]

896.3 [m]

896.3 [m]

896.3 [m]

1

[[XY]]

In most cases (14/18), the median daily distance of males from the nest they sired was <200 metres. There was one male whose median daily distance was slightly further (208 m), but he spent a large proportion of his time within 200 m (mean daily distance = 135 m).

The males with a greater median distance from the nest they sired have limited data (only one day).

We can also summarise the daily min, mean, median and max distance of successful males from nests they sired.

nest_male_date %>%
  filter(success==1) %>%
  ungroup() %>%
  summarise(min_dist=min(dist_m),
            max_dist=max(dist_m),
            mean_dist=mean(dist_m),
            median_dist=median(dist_m),
            sd_dist=sd(dist_m)) %>%
  flextable()

min_dist

max_dist

mean_dist

median_dist

sd_dist

geometry

8.7 [m]

1,350 [m]

210 [m]

116 [m]

270

[[XY]]

In a preliminary analysis (not included here), I explored the three cases where males sired nests far from their territory. P1103 was sired by a male who had a territory further away, but did visit the nest site.

P701 was sired by two males who had neighbouring territories on the other side of the site. The mother was seen on their territories when she was fertile, but apparently decided to build her nest somewhere else.

Determine likely paternity of predated nests

Males with a median daily distance <=200 metres from a nest site still seem to have a higher probability of siring the eggs. We use this information to get a list of “likely” fathers of the predated nests.

potential_fathers <- nest_male_date %>%
  filter(nest %in% predated_nests$nest) %>%
  group_by(nest, colour_combo) %>%
  summarise(first_dist = as.numeric(first(dist_m)),
            median_dist = as.numeric(median(dist_m)),
            mean_dist = as.numeric(mean(dist_m))) %>%
  filter(median_dist<=200) %>%
  ungroup() %>%
  group_by(colour_combo) %>%
  summarise(potential_n_extra_nests=n())

# to get rid of geometry from object
potential_fathers <- data.frame(colour_combo = potential_fathers$colour_combo,
                                potential_n_extra_nests = potential_fathers$potential_n_extra_nests)

potential_fathers %>% flextable()

colour_combo

potential_n_extra_nests

DB,R,O

2

DB,R,R

1

DG,PI,Y

1

DG,R,DG

1

O,DB,O

1

O,PI,Y

1

O,Y,DG

4

O,Y,O

2

PI,DG,DG

2

PI,O,O

2

PI,PI,PI

4

PI,R,DB

2

PI,Y,R

3

R,DG,R

1

R,O,R

1

R,Y,Y

2

Y,Y,Y

1

write.csv(potential_fathers, file=here("outputs","processed-data","potential_fathers.csv"), row.names=FALSE)
potential_fathers <- read.csv(here("outputs","processed-data","potential_fathers.csv"))

father_summary <- full_join(fathers,
                             potential_fathers,
                             by=join_by('LL.father'=='colour_combo')) %>%
  left_join(.,
            tagged_males[,c('colour_combo','bird_ID')],
            by=join_by('LL.father'=='colour_combo')) %>%
  ungroup() %>%
  select(!ID.father) %>%
  rename(ID.father=bird_ID) %>%
  mutate(success=case_when(n_eggs>0 ~ "yes",
                           TRUE ~ "no"),
         potential_n_extra_nests=replace_na(potential_n_extra_nests,0),
         likely_extra_nests=case_when(potential_n_extra_nests>0 ~ "yes",
                                       TRUE ~ "less likely"),
         likely_undetected_success=case_when(success=="no" & likely_extra_nests=="yes" ~ "yes",
                                             success=="yes" ~ "no",
                                        TRUE ~ "no")) %>%
  filter(!is.na(LL.father))

Combine success and activity data

activity <- left_join(activity,
                      father_summary,
                      by=join_by('bird_ID'=='ID.father',
                                 'colour_combo'=='LL.father')) %>%
  left_join(., tagged_males[,c('tag_ID','body_mass')], join_by('tag_ID')) %>%
  mutate(n_eggs=replace_na(n_eggs, 0),
         n_nests=replace_na(n_nests, 0),
         potential_n_extra_nests=replace_na(potential_n_extra_nests,0),
         likely_extra_nests=case_when(potential_n_extra_nests>0 ~ "yes",
                                      TRUE ~ "less likely"),
         success=case_when(n_eggs > 0 ~ "yes",
                           TRUE ~ "no"),
         likely_undetected_success=case_when(success=="no" & likely_extra_nests=="yes" ~ "yes",
                                        TRUE ~ "no"))

daily_activity <- left_join(daily_activity,
                      father_summary,
                      by=join_by('bird_ID'=='ID.father',
                                 'colour_combo'=='LL.father')) %>% 
    left_join(., tagged_males[,c('tag_ID','body_mass')], join_by('tag_ID')) %>%
  mutate(n_eggs=replace_na(n_eggs,0),
         n_nests=replace_na(n_nests, 0),
         potential_n_extra_nests=replace_na(potential_n_extra_nests,0),
         likely_extra_nests=case_when(potential_n_extra_nests>0 ~ "yes",
                                      TRUE ~ "less likely"),
         success=case_when(n_eggs > 0 ~ "yes",
                           TRUE ~ "no"),
         likely_undetected_success=case_when(success=="no" & likely_extra_nests=="yes" ~ "yes",
                                        TRUE ~ "no"))
activity <- activity %>%
  mutate(breeding_stage=case_when(window_end_localtime-10*60<=start_FertilePeriod ~ "pre-fertile",
                           window_end_localtime-10*60>start_FertilePeriod &
                             window_end_localtime<end_FertilePeriod ~ "fertile",
                           window_end_localtime>=end_FertilePeriod ~ "post-fertile"))

daily_activity <- daily_activity %>%
  mutate(breeding_stage=case_when(local_date<=start_FertilePeriod ~ "pre-fertile",
                           local_date>start_FertilePeriod &
                             local_date<end_FertilePeriod ~ "fertile",
                           local_date>=end_FertilePeriod ~ "post-fertile"))

daily_activity <- left_join(daily_activity, 
                             ff_by_date, 
                             by=join_by('local_date'=='date'),
                             relationship = "many-to-one") %>% 
  mutate(n_FertileFemales=replace_na(n_FertileFemales,0))

Summarise individual data and exclude incomplete data

daily_activity <- daily_activity %>%
  filter(n_hours>=min_acc_hrs_per_day)
ind_total_active <- daily_activity %>%
  group_by(bird_ID, colour_combo, breeding_stage,
           body_mass, n_eggs, n_nests, success,
           likely_undetected_success, potential_n_extra_nests, likely_extra_nests) %>%
  summarise(n_windows=sum(n_windows),
            n_days=n(),
            total_perc_active=sum(total_active)/sum(n_windows)*100)

# Note: total perc active will be same whether calculated from activity or daily_activity
# However, mean ODBA should be calculated from activity

ind_odba <- activity %>%
  filter(!is.na(activity_threshold)) %>%
    group_by(bird_ID, colour_combo, breeding_stage,
           body_mass, n_eggs, n_nests, success,
           likely_undetected_success, potential_n_extra_nests, likely_extra_nests) %>%
  summarise(mean_ODBA=mean(mean_ODBA))

ind_activity <- full_join(ind_total_active, ind_odba, 
                           by=c('bird_ID','colour_combo','breeding_stage',
                           'body_mass','n_eggs','n_nests','success',
                           'likely_undetected_success','potential_n_extra_nests',
                           'likely_extra_nests'))

# Filter out missing/insufficient data
ind_activity <- ind_activity %>% filter(!is.na(total_perc_active))


rm(ind_total_active, ind_odba)

Note: There are no males with <2 days of data during the fertile period and >=2 days of data during post-fertile period. None of the males with <2 days of data during the fertile period were successful.

males_with_sufficient_data = ind_activity %>%
  filter(breeding_stage=="fertile" & n_days>=min_acc_days) %>%
  select(colour_combo, bird_ID, n_days, n_windows, success)

males_with_sufficient_data %>%
  arrange(n_days, n_windows, success) %>%
  select(colour_combo, bird_ID, n_days, n_windows, n_nests) %>%
  flextable()

breeding_stage

body_mass

n_eggs

success

likely_undetected_success

potential_n_extra_nests

colour_combo

bird_ID

n_days

n_windows

n_nests

fertile

99.5

0

no

no

0

DB,PI,DB

141,269,115

2

268

0

fertile

97.4

0

no

no

0

Y,PI,DB

141,269,180

2

275

0

fertile

95.4

0

no

no

0

PI,O,PI

141,269,158

2

280

0

fertile

105.6

0

no

no

0

DG,PI,DB

141,269,037

2

288

0

fertile

105.7

0

no

no

0

R,R,Y

141,269,132

2

288

0

fertile

106.7

0

no

no

0

R,DB,O

141,269,039

3

396

0

fertile

90.4

0

no

no

0

R,Y,DB

141,269,159

3

417

0

fertile

106.4

0

no

no

0

DB,DB,O

141,269,031

3

431

0

fertile

109.8

0

no

no

0

O,R,Y

141,269,128

3

432

0

fertile

101.4

0

no

yes

2

R,Y,Y

141,269,173

3

432

0

fertile

96.6

4

yes

no

0

DG,DB,DG

141,269,127

3

432

1

fertile

93.7

0

no

no

0

PI,Y,PI

141,269,145

4

553

0

fertile

93.5

0

no

no

0

DG,Y,Y

141,269,038

4

569

0

fertile

96.7

0

no

no

0

DB,DB,Y

141,269,179

4

570

0

fertile

101.5

0

no

no

0

O,DG,R

141,269,004

4

576

0

fertile

99.3

0

no

no

0

DG,DB,PI

141,269,163

4

576

0

fertile

89.0

0

no

no

0

DG,PI,O

141,269,050

5

712

0

fertile

108.2

0

no

no

0

R,DG,Y

141,269,048

5

718

0

fertile

108.3

0

no

no

0

O,R,PI

141,269,047

5

720

0

fertile

105.5

0

no

no

0

Y,R,PI

141,269,151

5

720

0

fertile

103.5

0

no

yes

4

O,Y,DG

141,269,155

5

720

0

fertile

99.2

0

no

no

0

PI,DB,DB

141,269,183

5

720

0

fertile

95.3

4

yes

no

0

DB,DG,Y

141,269,046

5

720

1

fertile

95.2

0

no

no

0

R,O,O

141,269,045

6

843

0

fertile

97.9

2

yes

no

2

O,Y,O

141,269,182

6

845

1

fertile

92.5

5

yes

no

0

Y,DG,PI

141,269,044

6

846

2

fertile

93.2

0

no

no

0

R,Y,R

141,269,181

6

848

0

fertile

101.1

4

yes

no

0

DG,O,O

141,269,043

6

857

1

fertile

103.5

0

no

no

0

PI,PI,DB

141,269,178

6

861

0

fertile

105.7

0

no

no

0

R,Y,O

141,269,042

6

863

0

fertile

97.4

0

no

no

0

Y,Y,O

141,269,041

6

864

0

fertile

102.5

0

no

no

0

Y,DG,DG

141,269,167

6

864

0

fertile

97.6

0

no

no

0

DB,PI,PI

141,269,171

6

864

0

fertile

111.4

0

no

yes

1

DG,PI,Y

141,269,172

6

864

0

fertile

87.7

0

no

yes

4

PI,PI,PI

141,269,174

6

864

0

fertile

105.4

0

no

no

0

DB,R,Y

141,269,177

6

864

0

fertile

98.5

3

yes

no

0

DG,R,R

141,269,175

6

864

1

fertile

99.9

0

no

no

0

Y,Y,PI

141,269,040

7

987

0

fertile

108.7

0

no

no

0

R,O,DB

141,269,136

7

987

0

fertile

96.7

0

no

no

0

R,DB,R

141,269,170

7

989

0

fertile

97.5

4

yes

no

1

DB,R,R

141,269,169

7

1,000

1

fertile

97.9

0

no

no

0

DB,DG,DG

141,269,165

7

1,007

0

fertile

93.8

0

no

no

0

DB,DB,DB

141,269,006

7

1,008

0

fertile

95.7

0

no

no

0

DG,R,DB

141,269,034

7

1,008

0

fertile

102.5

0

no

no

0

R,PI,PI

141,269,036

7

1,008

0

fertile

102.4

0

no

no

0

DG,DG,DB

141,269,133

7

1,008

0

fertile

103.5

0

no

no

0

O,DB,R

141,269,160

7

1,008

0

fertile

115.9

0

no

no

0

DG,DB,Y

141,269,162

7

1,008

0

fertile

94.3

0

no

no

0

O,PI,O

141,269,164

7

1,008

0

fertile

100.4

0

no

yes

2

PI,O,O

141,269,166

7

1,008

0

fertile

107.3

0

no

no

0

DG,DG,R

141,269,033

8

1,144

0

fertile

104.7

0

no

yes

2

PI,R,DB

141,269,157

8

1,146

0

fertile

93.2

0

no

yes

1

O,PI,Y

141,269,156

8

1,150

0

fertile

103.9

0

no

no

0

DG,O,DG

141,269,030

8

1,152

0

fertile

105.5

0

no

yes

1

DG,R,DG

141,269,149

8

1,152

0

fertile

99.5

0

no

no

0

DG,R,O

141,269,153

8

1,152

0

fertile

100.4

0

no

no

0

O,R,O

141,269,154

8

1,152

0

fertile

108.1

3

yes

no

2

DB,R,O

141,269,148

8

1,152

1

fertile

100.5

1

yes

no

1

Y,Y,Y

141,269,152

8

1,152

1

fertile

105.2

5

yes

no

0

O,DG,DG

141,269,147

9

1,276

2

fertile

101.8

2

yes

no

0

R,DG,PI

141,269,146

9

1,281

1

fertile

102.4

3

yes

no

0

DB,Y,R

141,269,144

9

1,289

1

fertile

102.2

0

no

no

0

Y,O,R

141,269,143

9

1,292

0

fertile

106.9

0

no

yes

1

R,DG,R

141,269,137

9

1,295

0

fertile

97.5

0

no

yes

1

R,O,R

141,269,138

9

1,296

0

fertile

89.2

0

no

no

0

PI,O,R

141,269,140

9

1,296

0

fertile

88.4

0

no

no

0

R,R,DG

141,269,141

9

1,296

0

fertile

94.7

4

yes

no

1

O,DB,O

141,269,026

9

1,296

1

fertile

97.2

0

no

yes

3

PI,Y,R

141,269,011

10

1,422

0

fertile

99.6

0

no

no

0

DB,Y,PI

141,269,134

10

1,436

0

fertile

109.6

0

no

no

0

PI,PI,O

141,269,024

10

1,437

0

fertile

97.5

0

no

no

0

Y,DB,Y

141,269,022

10

1,440

0

fertile

104.3

0

no

no

0

Y,Y,DG

141,269,129

10

1,440

0

fertile

103.2

0

no

no

0

Y,O,Y

141,269,126

11

1,566

0

fertile

99.2

0

no

yes

2

PI,DG,DG

141,269,122

11

1,584

0

fertile

95.8

1

yes

no

0

Y,Y,R

141,269,125

11

1,584

1

fertile

93.5

0

no

no

0

PI,Y,Y

141,269,105

12

1,720

0

fertile

100.1

0

no

no

0

O,DB,DG

141,269,010

12

1,728

0

fertile

107.5

0

no

no

0

DB,O,DG

141,269,118

12

1,728

0

fertile

101.4

6

yes

no

0

PI,PI,R

141,269,009

12

1,728

2

fertile

106.0

3

yes

no

0

PI,R,PI

141,269,114

13

1,871

1

There are 81 males with at least 2 days of complete (20.4 hours) data during the fertile period, as shown in the table above.

From here on, the activity datasets are restricted only to these males.

daily_activity <- daily_activity %>%
  filter(bird_ID %in% males_with_sufficient_data$bird_ID &
           n_hours>=min_acc_hrs_per_day)

ind_activity <- ind_activity %>%
  filter(bird_ID %in% males_with_sufficient_data$bird_ID)
ind_activity <- left_join(ind_activity,
                          male_tenure,
                          join_by('bird_ID', 'colour_combo'))

Q: Do males spend more time active when there are more mating opportunities?

The percentage of time that males spent active in a day ranged from less than 20% to over 99% (mean ± standard deviation: 55.634 +/- 21.494.

daily_activity %>%
  ungroup() %>%
  summarise(min_perc_active=min(perc_active),
            max_perc_active=max(perc_active),
            mean_perc_active=mean(perc_active),
            sd_perc_active=sd(perc_active),
            n_birds=n_distinct(bird_ID),
            n=n(),
            se_perc_active=sd_perc_active/sqrt(n)) %>%
  flextable()

min_perc_active

max_perc_active

mean_perc_active

sd_perc_active

n_birds

n

se_perc_active

4.86

100

55.6

21.5

81

802

0.759

within_ind_var <- daily_activity %>%
  ungroup() %>%
  select(colour_combo, perc_active, mean_ODBA) %>%
  group_by(colour_combo) %>%
  summarise(min_perc_active=min(perc_active),
            max_perc_active=max(perc_active),
            diff_perc_active=max_perc_active-min_perc_active,
            diff_hours_active=diff_perc_active*24/100,
            diff_multiply_active=max_perc_active/min_perc_active,
            mean_perc_active=mean(perc_active),
            sd_perc_active=sd(perc_active),
            min_ODBA=min(mean_ODBA),
            max_ODBA=max(mean_ODBA),
            diff_ODBA=max_ODBA-min_ODBA,
            diff_multiply_ODBA=max_ODBA/min_ODBA,
            mean_mean_ODBA=mean(mean_ODBA),
            sd_ODBA=sd(mean_ODBA),
            n_days=n()) %>%
  arrange(diff_multiply_ODBA)

within_ind_var %>% flextable()

colour_combo

min_perc_active

max_perc_active

diff_perc_active

diff_hours_active

diff_multiply_active

mean_perc_active

sd_perc_active

min_ODBA

max_ODBA

diff_ODBA

diff_multiply_ODBA

mean_mean_ODBA

sd_ODBA

n_days

Y,PI,DB

55.73

63.2

7.47

1.793

1.13

59.5

5.28

0.2369

0.250

0.0129

1.05

0.243

0.00913

2

DB,PI,DB

58.06

60.4

2.35

0.565

1.04

59.2

1.66

0.3175

0.347

0.0292

1.09

0.332

0.02061

2

O,DG,R

63.19

75.0

11.81

2.833

1.19

69.6

5.93

0.3069

0.345

0.0377

1.12

0.333

0.01768

4

DB,DB,DB

49.31

63.9

14.58

3.500

1.30

58.1

6.12

0.2453

0.308

0.0631

1.26

0.283

0.02566

7

PI,Y,PI

48.09

77.1

28.99

6.958

1.60

60.5

12.14

0.2521

0.323

0.0708

1.28

0.280

0.03024

4

O,PI,O

39.58

54.9

15.28

3.667

1.39

48.5

5.68

0.1944

0.251

0.0564

1.29

0.222

0.02284

7

R,Y,DB

31.25

63.9

32.64

7.833

2.04

48.3

16.36

0.2769

0.359

0.0820

1.30

0.305

0.04651

3

DG,PI,DB

38.19

54.2

15.97

3.833

1.42

46.2

11.29

0.1918

0.249

0.0576

1.30

0.221

0.04075

2

Y,R,PI

41.67

62.5

20.83

5.000

1.50

53.2

7.63

0.1424

0.186

0.0436

1.31

0.165

0.01587

5

Y,O,Y

27.08

52.8

25.69

6.167

1.95

41.1

6.96

0.1656

0.217

0.0517

1.31

0.189

0.01862

13

Y,Y,O

40.28

61.1

20.83

5.000

1.52

51.5

6.80

0.1654

0.223

0.0572

1.35

0.198

0.01898

10

DG,DB,DG

36.11

59.0

22.92

5.500

1.63

49.3

11.85

0.1972

0.268

0.0708

1.36

0.240

0.03756

3

DB,DB,O

46.53

67.1

20.61

4.945

1.44

57.8

10.44

0.2191

0.303

0.0840

1.38

0.256

0.04302

3

DB,Y,PI

29.86

47.9

18.06

4.333

1.60

39.4

6.42

0.1486

0.212

0.0633

1.43

0.178

0.02132

12

DG,R,DB

40.97

66.7

25.69

6.167

1.63

56.2

6.98

0.1841

0.268

0.0839

1.46

0.231

0.02551

9

R,DG,Y

36.62

59.0

22.41

5.378

1.61

49.6

8.20

0.1663

0.243

0.0763

1.46

0.210

0.02835

10

O,R,Y

52.78

64.6

11.81

2.833

1.22

60.4

6.62

0.2075

0.303

0.0954

1.46

0.257

0.04783

3

R,R,Y

27.78

51.4

23.61

5.667

1.85

39.6

16.70

0.1542

0.232

0.0782

1.51

0.193

0.05529

2

DB,DB,Y

58.33

91.3

32.97

7.913

1.57

73.7

14.02

0.3353

0.512

0.1766

1.53

0.417

0.07697

4

DG,DB,Y

38.19

65.3

27.08

6.500

1.71

49.7

10.69

0.1642

0.253

0.0892

1.54

0.200

0.03196

9

Y,DG,DG

31.25

56.9

25.69

6.167

1.82

47.3

8.72

0.1509

0.234

0.0835

1.55

0.207

0.02755

7

DG,DB,PI

39.58

70.1

30.56

7.333

1.77

56.4

13.71

0.1915

0.300

0.1086

1.57

0.250

0.05316

4

DG,R,O

26.39

50.0

23.61

5.667

1.89

39.1

8.33

0.1194

0.188

0.0681

1.57

0.146

0.02204

11

DB,R,Y

23.61

54.2

30.56

7.333

2.29

37.3

9.10

0.1300

0.206

0.0758

1.58

0.162

0.02303

9

PI,O,PI

4.86

26.5

21.61

5.186

5.45

15.7

15.28

0.1356

0.217

0.0812

1.60

0.176

0.05741

2

DG,PI,O

31.94

53.4

21.44

5.145

1.67

43.9

8.02

0.1389

0.224

0.0851

1.61

0.177

0.02543

9

PI,PI,DB

14.18

39.6

25.40

6.096

2.79

34.5

10.05

0.1099

0.179

0.0688

1.63

0.160

0.02548

6

R,Y,Y

39.58

73.6

34.03

8.167

1.86

52.1

11.31

0.1779

0.294

0.1157

1.65

0.225

0.04828

7

R,DB,O

17.83

49.6

31.76

7.623

2.78

31.0

16.54

0.1331

0.221

0.0874

1.66

0.174

0.04393

3

R,O,DB

19.51

48.6

29.10

6.984

2.49

39.4

10.29

0.1099

0.183

0.0733

1.67

0.162

0.02591

7

DG,Y,Y

42.66

75.0

32.34

7.762

1.76

57.8

13.43

0.2381

0.400

0.1614

1.68

0.314

0.06655

4

DB,PI,PI

40.28

75.0

34.72

8.333

1.86

55.6

11.11

0.1386

0.241

0.1020

1.74

0.181

0.02787

14

O,PI,Y

22.92

58.3

35.42

8.500

2.55

39.9

10.29

0.1375

0.239

0.1017

1.74

0.181

0.02529

24

PI,DB,DB

30.37

75.7

45.32

10.878

2.49

50.6

15.24

0.1575

0.275

0.1173

1.74

0.204

0.03921

10

O,R,PI

44.44

77.8

33.33

8.000

1.75

62.4

12.27

0.1922

0.336

0.1437

1.75

0.260

0.04871

8

DG,R,R

31.94

72.2

40.28

9.667

2.26

52.2

13.64

0.2296

0.402

0.1727

1.75

0.313

0.06048

8

PI,O,R

36.11

72.9

36.81

8.833

2.02

50.6

12.38

0.1767

0.328

0.1514

1.86

0.228

0.05180

11

PI,Y,R

11.11

59.0

47.92

11.500

5.31

37.7

11.51

0.1314

0.245

0.1140

1.87

0.193

0.02576

20

Y,DG,PI

29.17

68.8

39.58

9.500

2.36

48.5

11.69

0.1375

0.258

0.1209

1.88

0.193

0.03527

16

DG,DG,DB

52.78

98.6

45.83

11.000

1.87

78.0

19.35

0.3233

0.682

0.3584

2.11

0.498

0.15002

7

DG,R,DG

19.44

68.8

49.31

11.833

3.54

51.0

13.63

0.1174

0.254

0.1365

2.16

0.202

0.03907

17

O,Y,DG

66.67

99.3

32.64

7.833

1.49

81.5

12.41

0.3059

0.663

0.3574

2.17

0.432

0.14224

5

Y,Y,Y

58.33

99.3

40.97

9.833

1.70

79.4

15.48

0.2734

0.602

0.3288

2.20

0.409

0.12318

10

O,DB,O

37.50

91.7

54.17

13.000

2.44

63.4

18.19

0.2571

0.574

0.3171

2.23

0.381

0.09446

11

R,R,DG

50.69

100.0

49.31

11.833

1.97

73.7

15.03

0.2520

0.583

0.3306

2.31

0.371

0.09908

11

PI,O,O

36.81

97.9

61.11

14.667

2.66

67.4

21.93

0.2311

0.566

0.3351

2.45

0.353

0.12218

7

DB,Y,R

31.94

91.7

59.72

14.333

2.87

57.1

18.44

0.1691

0.414

0.2454

2.45

0.258

0.08187

12

PI,R,DB

46.53

100.0

53.47

12.833

2.15

73.7

17.86

0.2597

0.642

0.3826

2.47

0.409

0.13409

9

R,DB,R

34.03

84.0

50.00

12.000

2.47

57.3

12.46

0.1640

0.413

0.2491

2.52

0.252

0.06536

15

DG,O,O

35.42

83.9

48.52

11.646

2.37

59.4

14.54

0.1905

0.500

0.3092

2.62

0.294

0.08758

10

DB,DG,DG

20.14

90.2

70.07

16.817

4.48

49.9

20.88

0.1682

0.441

0.2733

2.62

0.282

0.08086

10

R,Y,O

33.33

93.0

59.67

14.322

2.79

55.8

23.18

0.1753

0.462

0.2869

2.64

0.270

0.11065

6

R,PI,PI

32.64

96.5

63.89

15.333

2.96

63.0

22.96

0.2363

0.634

0.3977

2.68

0.376

0.13598

8

PI,R,PI

18.75

89.6

70.83

17.000

4.78

66.0

18.05

0.2153

0.580

0.3645

2.69

0.400

0.08569

16

PI,PI,O

36.11

98.6

62.50

15.000

2.73

70.0

18.53

0.2158

0.581

0.3654

2.69

0.369

0.10271

14

O,DB,R

9.03

71.5

62.50

15.000

7.92

47.8

17.05

0.1008

0.271

0.1707

2.69

0.197

0.05015

11

Y,Y,R

22.92

94.4

71.53

17.167

4.12

61.9

25.87

0.1849

0.498

0.3131

2.69

0.320

0.12107

13

R,DG,PI

39.58

98.6

59.03

14.167

2.49

70.1

17.17

0.2219

0.599

0.3774

2.70

0.368

0.11371

10

PI,PI,R

5.56

100.0

94.44

22.667

18.00

50.4

31.79

0.2097

0.604

0.3944

2.88

0.391

0.13949

15

DG,O,DG

30.56

97.2

66.67

16.000

3.18

63.8

23.44

0.1700

0.495

0.3248

2.91

0.301

0.10971

11

Y,O,R

15.97

92.4

76.39

18.333

5.78

64.4

22.35

0.1753

0.513

0.3372

2.92

0.363

0.09593

14

PI,Y,Y

31.62

98.6

66.99

16.078

3.12

76.1

19.28

0.2359

0.699

0.4632

2.96

0.489

0.13686

12

DB,O,DG

7.64

95.8

88.19

21.167

12.55

44.7

25.67

0.1652

0.506

0.3409

3.06

0.293

0.10755

22

DB,DG,Y

14.58

79.2

64.58

15.500

5.43

35.8

20.08

0.1253

0.390

0.2649

3.11

0.191

0.08119

10

O,DB,DG

17.36

98.6

81.25

19.500

5.68

62.6

26.71

0.1691

0.536

0.3670

3.17

0.348

0.12297

22

Y,DB,Y

18.75

92.4

73.61

17.667

4.93

53.5

25.28

0.1608

0.516

0.3553

3.21

0.284

0.12621

18

R,O,O

22.22

88.6

66.40

15.935

3.99

43.0

22.17

0.1052

0.341

0.2355

3.24

0.178

0.07560

9

O,Y,O

25.69

88.0

62.31

14.953

3.42

61.7

24.14

0.1528

0.498

0.3455

3.26

0.334

0.13877

8

O,DG,DG

25.00

98.6

73.60

17.664

3.94

73.9

21.22

0.1613

0.528

0.3664

3.27

0.364

0.10335

15

DB,R,R

38.89

97.8

58.91

14.137

2.51

67.8

20.09

0.1943

0.642

0.4478

3.30

0.341

0.12866

11

PI,DG,DG

31.94

98.6

66.67

16.000

3.09

64.5

24.22

0.1797

0.599

0.4192

3.33

0.350

0.15376

16

DG,PI,Y

36.11

98.6

62.50

15.000

2.73

61.5

21.20

0.1895

0.637

0.4472

3.36

0.323

0.14337

10

PI,PI,PI

24.31

96.5

72.22

17.333

3.97

67.2

19.11

0.1601

0.552

0.3918

3.45

0.339

0.09913

11

O,R,O

34.72

98.6

63.89

15.333

2.84

65.3

23.69

0.1588

0.548

0.3891

3.45

0.309

0.14162

10

DG,DG,R

24.31

97.9

73.61

17.667

4.03

58.4

26.36

0.1588

0.557

0.3979

3.51

0.305

0.13729

12

R,O,R

37.50

97.9

60.42

14.500

2.61

66.7

21.59

0.1706

0.601

0.4301

3.52

0.334

0.14649

12

Y,Y,PI

25.00

96.0

70.97

17.032

3.84

50.8

25.68

0.1199

0.428

0.3077

3.57

0.225

0.12562

8

Y,Y,DG

15.97

97.9

81.94

19.667

6.13

44.3

27.91

0.1535

0.655

0.5020

4.27

0.277

0.15214

20

DB,R,O

11.81

99.3

87.50

21.000

8.41

61.6

29.07

0.1249

0.578

0.4536

4.63

0.316

0.15355

14

R,Y,R

14.58

90.3

75.69

18.167

6.19

49.2

28.02

0.0967

0.505

0.4078

5.22

0.245

0.14058

14

R,DG,R

19.44

100.0

80.56

19.333

5.14

60.0

31.50

0.1416

0.807

0.6657

5.70

0.316

0.20293

12

daily_activity %>%
  ungroup() %>%
  group_by(breeding_stage) %>%
  summarise(mean_perc_active=mean(perc_active),
            sd_perc_active=sd(perc_active),
            n_samples=n(),
            n_birds=n_distinct(bird_ID),
            se_perc_active=sd_perc_active/sqrt(n_samples)) %>%
  flextable()

breeding_stage

mean_perc_active

sd_perc_active

n_samples

n_birds

se_perc_active

fertile

62.1

20.7

555

81

0.877

post-fertile

41.0

15.3

247

56

0.972

daily_activity %>%
  ungroup() %>%
  group_by(breeding_stage) %>%
  summarise(min_mean_ODBA=min(mean_ODBA),
            max_mean_ODBA=max(mean_ODBA),
            mean_mean_ODBA=mean(mean_ODBA),
            sd_mean_ODBA=sd(mean_ODBA),
            n_samples=n(),
            n_birds=n_distinct(bird_ID),
            se_mean_ODBA=sd_mean_ODBA/sqrt(n_samples)) %>%
  flextable()

breeding_stage

min_mean_ODBA

max_mean_ODBA

mean_mean_ODBA

sd_mean_ODBA

n_samples

n_birds

se_mean_ODBA

fertile

0.1008

0.807

0.315

0.1320

555

81

0.0056

post-fertile

0.0967

0.467

0.205

0.0518

247

56

0.0033

daily_activity %>%
  ungroup() %>%
  select(colour_combo, perc_active, mean_ODBA, breeding_stage) %>%
  group_by(colour_combo, breeding_stage) %>%
  summarise(mean_perc_active=mean(perc_active),
            mean_ODBA=mean(mean_ODBA)) %>%
  pivot_wider(values_from = c('mean_perc_active','mean_ODBA'),
              names_from = 'breeding_stage') %>%
  filter(!is.na(`mean_ODBA_post-fertile`) & !is.na(mean_perc_active_fertile)) %>%
  mutate(diff_perc_active=mean_perc_active_fertile-`mean_perc_active_post-fertile`,
         diff_hours_active=24/100*diff_perc_active,
         diff_mean_ODBA=mean_ODBA_fertile-`mean_ODBA_post-fertile`) %>%
  ungroup() %>%
  summarise(mean_diff_perc_active=mean(diff_perc_active),
            sd_diff_perc_active=sd(diff_perc_active),
            mean_diff_hours_active=mean(diff_hours_active),
            sd_diff_hours_active=sd(diff_hours_active),
            mean_diff_ODBA=mean(diff_mean_ODBA),
            sd_diff_ODBA=sd(diff_mean_ODBA)) %>%
  flextable()

mean_diff_perc_active

sd_diff_perc_active

mean_diff_hours_active

sd_diff_hours_active

mean_diff_ODBA

sd_diff_ODBA

20.6

18.1

4.95

4.35

0.0991

0.0832

daily_activity %>%
  ungroup() %>%
  filter(n_FertileFemales==0 | n_FertileFemales==40) %>%
  select(colour_combo, perc_active, n_FertileFemales) %>%
  group_by(colour_combo, n_FertileFemales) %>%
  summarise(mean_perc_active=mean(perc_active)) %>%
  pivot_wider(values_from = 'mean_perc_active', names_from = 'n_FertileFemales', names_prefix="ff_") %>%
  filter(!is.na(ff_40) & !is.na(ff_0)) %>%
  mutate(diff=ff_40-ff_0,
         diff_hours=24/100*diff) %>%
  flextable()

colour_combo

ff_40

ff_0

diff

diff_hours

DB,O,DG

81.2

20.0

61.23

14.694

DB,PI,PI

54.9

62.8

-7.99

-1.917

DB,R,O

99.0

19.4

79.51

19.083

DG,R,DG

53.8

49.9

3.90

0.936

O,DB,DG

95.5

28.8

66.67

16.000

O,DG,DG

97.9

28.8

69.09

16.582

O,PI,Y

35.9

39.1

-3.20

-0.768

PI,DG,DG

98.6

34.0

64.58

15.500

PI,PI,PI

96.5

24.3

72.22

17.333

PI,Y,R

38.2

25.6

12.62

3.028

R,DB,R

78.0

56.9

21.07

5.057

R,Y,R

66.4

22.1

44.27

10.624

Y,DB,Y

85.1

30.0

55.03

13.208

Y,DG,PI

57.9

46.1

11.87

2.849

Y,O,R

89.6

16.0

73.61

17.667

Y,Y,DG

69.1

22.9

46.18

11.083

The table above shows the percentage of time that males (each with a unique colour_combo) were active when female fertility was at its peak (40 fertile females, ‘ff_40’) compared with when females were no longer fertile and were incubating (0 fertile females, ‘ff_0’). Many males were much more active during the peak fertile period.

The table below summarises activity levels when there were no fertile females.

daily_activity %>%
  ungroup() %>%
  filter(n_FertileFemales==0) %>%
  summarise(mean_perc_active=mean(perc_active),
            sd_perc_active=sd(perc_active),
            mean_hours_active=24/100*mean_perc_active,
            sd_hours_active=24/100*sd_perc_active) %>%
  flextable()

mean_perc_active

sd_perc_active

mean_hours_active

sd_hours_active

34.7

15.6

8.33

3.74

The table below summarises activity levels when there were the maximum number of fertile females.

daily_activity %>%
  ungroup() %>%
  filter(n_FertileFemales==40) %>%
  summarise(mean_perc_active=mean(perc_active),
            sd_perc_active=sd(perc_active),
            mean_hours_active=24/100*mean_perc_active,
            sd_hours_active=24/100*sd_perc_active) %>%
  flextable()

mean_perc_active

sd_perc_active

mean_hours_active

sd_hours_active

71.8

25.4

17.2

6.1

daily_activity %>%
  ungroup() %>%
  filter(n_FertileFemales==0 | n_FertileFemales==40) %>%
  select(colour_combo, perc_active, n_FertileFemales) %>%
  group_by(colour_combo, n_FertileFemales) %>%
  summarise(mean_perc_active=mean(perc_active)) %>%
  pivot_wider(values_from = 'mean_perc_active', names_from = 'n_FertileFemales', names_prefix="ff_") %>%
  filter(!is.na(ff_40) & !is.na(ff_0)) %>%
  mutate(diff=ff_40-ff_0,
         diff_hours=24/100*diff) %>%
  ungroup() %>%
  summarise(average_within_ind_diff_hours = mean(diff_hours),
            sd_within_ind_diff_hours = sd(diff_hours),
            n=n()) %>%
  flextable()

average_within_ind_diff_hours

sd_within_ind_diff_hours

n

10.1

7.32

16

Within individual males, the average difference in activity between peak fertile females and no fertile females was over 10 hours.

daily_activity %>%
  ungroup() %>%
  filter(n_FertileFemales==0 | n_FertileFemales==40) %>%
  select(colour_combo, perc_active, n_FertileFemales) %>%
  group_by(n_FertileFemales) %>%
  summarise(mean_perc_active=mean(perc_active)) %>%
  pivot_wider(values_from = 'mean_perc_active', names_from = 'n_FertileFemales', names_prefix="ff_") %>%
  mutate(diff=ff_40-ff_0,
         diff_hours=24/100*diff) %>%
  ungroup() %>%
  summarise(average_among_ind_diff_hours = mean(diff_hours)) %>%
  flextable()

average_among_ind_diff_hours

8.91

Among all males, the average difference was just under 9 hours.

Prediction: Males increase their activity (% time active per day) when there are fertile females at the site.

postfertile_males <- ind_activity %>%
  filter(breeding_stage=="post-fertile" & n_days>=2) %>%
  group_by(bird_ID) %>%
  summarise(mean(total_perc_active),
            mean(mean_ODBA))

fertile_males <- ind_activity %>%
  filter(breeding_stage=="fertile") %>% # already filtered to n_days>=2
  group_by(bird_ID) %>%
  summarise(mean(total_perc_active),
            mean(mean_ODBA))

Fig: Breeding stage and male % active

ind_activity_sub <- ind_activity %>%
  filter(n_days>=2 & breeding_stage!="pre-fertile") %>%
  mutate(breeding_stage=as.factor(breeding_stage))

ind_activity_plotdat <- rbind(ind_activity_sub %>% mutate(plot_type="box"), 
                              ind_activity_sub %>% mutate(plot_type="point"))

ind_activity_plotdat$plot_factor <- factor(
  paste(ind_activity_plotdat$breeding_stage, ind_activity_plotdat$plot_type),
  levels = c("fertile box", "fertile point", "post-fertile point", "post-fertile box"))

activity_fert_boxplot <- ggplot() +
  theme_classic()+
  geom_boxplot(data = ind_activity_plotdat %>% 
                 filter(plot_factor == "fertile box"),
               aes(x = plot_factor,
                   y = total_perc_active),
               outlier.shape=default_shape,
               fill="white",
               width=0.7) +
  geom_point(data = ind_activity_plotdat %>% 
                 filter(plot_factor == "fertile point" |
                          plot_factor == "post-fertile point"),
               aes(x = plot_factor,
                   y = total_perc_active,
                   fill=success, colour=success, alpha=success),
             size=small.point.size,
             position=position_jitter(w=0.15, h=0, seed = 1)) +
  geom_boxplot(data = ind_activity_plotdat %>% 
                 filter(plot_factor == "post-fertile box"),
               aes(x = plot_factor,
                   y = total_perc_active),
               outlier.shape=default_shape,
               fill="white",
               width=0.7)+
  geom_line(data = ind_activity_plotdat %>% 
                 filter(plot_factor == "fertile point" |
                          plot_factor == "post-fertile point"),
            aes(x = plot_factor,
                y = total_perc_active,
                colour=success, group=bird_ID, alpha=success), 
            size=small.point.size/3,
            position=position_jitter(w=0.15, h=0, seed = 1)) +
  coord_cartesian(ylim=c(0,100))+
  scale_y_continuous(breaks=c(seq(0,100,by=20)),
                     expand=c(0,0))+
  scale_x_discrete(breaks = c("fertile box", "post-fertile box"),
                   labels = c("fertile","post-fertile"))+
  scale_colour_manual(values=c(fail_colour,success_colour), name="Sired offspring?")+
  scale_fill_manual(values=c(fail_colour,success_colour), name="Sired offspring?")+
  scale_alpha_manual(values=c(0.5,0.8), name="Sired offspring?")+
  scale_shape_manual(values=c(default_shape), guide="none") +
  xlab("\n breeding stage")+
  ylab("male activity (% time)\n")+
  theme(text=element_text(family=font.family),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size),
        legend.title = element_text(size=legend.heading.size),
        legend.text = element_text(size=label.size),
        legend.position="right",
        legend.justification="top",
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "grey"))

activity_fert_boxplot

ind_activity %>%
  group_by(breeding_stage, success) %>%
  summarise(mean_total_perc_active = mean(total_perc_active),
            sd_total_perc_active = sd(total_perc_active),
            n=n(),
            se_total_perc_active = sd_total_perc_active/sqrt(n)) %>%
  flextable()

breeding_stage

success

mean_total_perc_active

sd_total_perc_active

n

se_total_perc_active

fertile

no

58.2

14.9

65

1.85

fertile

yes

67.1

12.1

16

3.03

post-fertile

no

41.9

10.5

41

1.64

post-fertile

yes

40.3

12.6

15

3.25

In the plot above, there are 81 individual males in total, PI,PI,R, O,DB,DG, PI,Y,R, Y,DB,Y, PI,PI,O, O,DB,O, DG,O,DG, DG,DG,R, DG,R,DB, R,PI,PI, Y,Y,PI, Y,Y,O, DG,O,O, Y,DG,PI, R,O,O, DB,DG,Y, O,R,PI, R,DG,Y, DG,PI,O, PI,R,PI, DB,O,DG, PI,DG,DG, Y,Y,R, Y,O,Y, Y,Y,DG, DB,Y,PI, R,DG,R, R,O,R, PI,O,R, R,R,DG, Y,O,R, DB,Y,R, R,DG,PI, O,DG,DG, DB,R,O, DG,R,DG, Y,Y,Y, DG,R,O, O,R,O, O,PI,Y, PI,R,DB, O,DB,R, DG,DB,Y, DB,DG,DG, Y,DG,DG, DB,R,R, R,DB,R, DB,PI,PI, DG,PI,Y, R,Y,Y, PI,PI,PI, DG,R,R, DB,R,Y, R,Y,R, O,Y,O, PI,DB,DB of which have post-fertile data. Males with “uncertain” success are treated as “no” success. Their inclusion or exclusion does not ever seem to affect results.

Each datapoint represents the mean activity of an individual male. The pattern we see here is similar to what was found in Lesku et al. 2012, although activity levels are shifted around 30% lower.

Model: Breeding stage and male % active

daily_activity <- daily_activity %>%
  mutate(prop_activity=perc_active/100,
         bird_ID=as.factor(bird_ID),
         breeding_stage=as.factor(breeding_stage),
         success=as.factor(success)) %>%
  mutate(prop_activity=case_when(prop_activity==1 ~ 0.99999,
                               TRUE ~ prop_activity))
stage_m <- glmmTMB(prop_activity ~ breeding_stage + (breeding_stage|bird_ID),
                   data=daily_activity,
                   #data=daily_activity %>% filter(likely_undetected_success=="no"),
                   family=beta_family(link="logit"))
summary(stage_m)
##  Family: beta  ( logit )
## Formula:          prop_activity ~ breeding_stage + (breeding_stage | bird_ID)
## Data: daily_activity
## 
##      AIC      BIC   logLik deviance df.resid 
##     -518     -490      265     -530      796 
## 
## Random effects:
## 
## Conditional model:
##  Groups  Name                       Variance Std.Dev. Corr  
##  bird_ID (Intercept)                0.415    0.644          
##          breeding_stagepost-fertile 0.585    0.765    -0.93 
## Number of obs: 802, groups:  bird_ID, 81
## 
## Dispersion parameter for beta family (): 6.78 
## 
## Conditional model:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                  0.5391     0.0796    6.77  1.2e-11 ***
## breeding_stagepost-fertile  -0.8601     0.1073   -8.02  1.1e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(stage_m, level=0.95)
##                                                     2.5 % 97.5 % Estimate
## (Intercept)                                         0.383  0.695    0.539
## breeding_stagepost-fertile                         -1.070 -0.650   -0.860
## Std.Dev.(Intercept)|bird_ID                         0.532  0.780    0.644
## Std.Dev.breeding_stagepost-fertile|bird_ID          0.605  0.968    0.765
## Cor.breeding_stagepost-fertile.(Intercept)|bird_ID -0.967 -0.774   -0.929
simulationOutput <- simulateResiduals(stage_m)
testResiduals(simulationOutput)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.06, p-value = 0.004
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.8, p-value = 0.02
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 7, observations = 802, p-value = 0.7
## alternative hypothesis: true probability of success is not equal to 0.00797
## 95 percent confidence interval:
##  0.00352 0.01790
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                0.00873

Important note We are using the DHARMA package to visualise whether models met model assumptions. This package also provides statistical tests and p-values for assumption violations, but we have focused on the visualisation of these assumptions instead. Sometimes there were statistically significant deviations from a particular assumption, but these appear reasonably minor and we believe they are of no cause for concern.

stage_m_flex <- as_flextable(stage_m)
stage_m_flex

group

Estimate

Standard Error

statistic

p-value

Fixed effects

(Intercept)

0.539

0.080

6.775

0.0000

***

breeding_stagepost-fertile

-0.860

0.107

-8.016

0.0000

***

Random effects

bird_ID

sd__(Intercept)

0.644

bird_ID

sd__breeding_stagepost-fertile

0.765

bird_ID

cor__(Intercept).breeding_stagepost-fertile

-0.929

Signif. codes: 0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

square root of the estimated residual variance: 6.8

data's log-likelihood under the model: 265.2

Akaike Information Criterion: -518.5

Bayesian Information Criterion: -490.4

if(save_tables==TRUE){
save_as_docx(stage_m_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/breeding_stage_model.docx"))
}

When ‘potentially’ successful males are excluded the results remain the same. There are no significant issues with model dispersion or influence of outliers.

Male spent more time active during the fertile stage, compared with the post-fertile stage. The inclusion of males that were more likely to have sired additional offspring (predated before they could be sampled) had no effect on these results.

Prediction 2 Male activity (% time active per day) increases with the number of fertile females at the site.

We also examine this trend in more detail, by looking at activity and the number of fertile females each day across the season.

Fig: Date and male % active

date_activity_n <- daily_activity %>%
  group_by(local_date) %>%
  summarise(n_ind=n_distinct(colour_combo))
date_labels = c("",
                "12 June", "", "", "",
                "16 June", "", "", "",
                "20 June", "", "", "",
                "24 June", "", "", "",
                "28 June", "", "", "",
                "2 July", "","","",
                "6 July", "","","",
                "10 July", "")
boxplot_active_date <- ggplot(daily_activity, aes(x=local_date,
                                                 y=perc_active,
                                                 group=local_date))+
  # theme_classic()+
  geom_boxplot(colour="gray30", outlier.shape = NA)+
  geom_point(size=small.point.size,
             aes(colour=success, alpha=success, fill=success),
             shape=default_shape)+
  geom_vline(xintercept=c(seq(min(daily_activity$local_date)-24*60*30,
                              max(daily_activity$local_date)+24*60*30,
                              by="days")),
             size=0.5,
             colour="darkgrey")+
  geom_vline(xintercept=max(daily_activity$local_date)+24*60*30,
             size=0.5,
             colour="black")+
  coord_flip(ylim=c(0,115),
             xlim = c(max(daily_activity$local_date)+24*60*30,
                      min(daily_activity$local_date)-24*60*30))+
  scale_x_continuous(breaks = c(seq(
                              local_date_to_posix('2022-06-11'),
                              local_date_to_posix('2022-07-11'),
                              by="days")),
                     labels=date_labels,
                     expand=c(0,0))+
  scale_y_continuous(breaks=seq(0,100,by=20),
                     expand=c(0,0)) +
  scale_colour_manual(values = c(fail_colour, success_colour), guide="none")+
  scale_fill_manual(values = c(fail_colour, success_colour), guide="none")+
  scale_alpha_manual(values = c(0.5, 0.8), guide="none")+
  xlab("date\n")+
  ylab("\nmale activity (% time)")+
  geom_text(data = date_activity_n,
            family=font.family,
            size=(label.size-2)/.pt,
            colour="gray50",
            aes(local_date, 107, label = n_ind),
            vjust = 0.3,
            hjust = 0)+
  theme(axis.ticks.length.y = unit(0.07, "cm"),
        text=element_text(family=font.family),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size),
        axis.line.y = element_line(size=0.5, color="black"),
        panel.background=element_blank())

boxplot_active_date

hor_barplot_ff <- barplot_ff +
  theme(axis.line.x=element_blank()) +
  coord_flip(xlim = c(max(daily_activity$local_date)+24*60*30, 
                      min(daily_activity$local_date)-24*60*30),
             ylim=c(0,42)) +
  geom_vline(xintercept=c(seq(min(daily_activity$local_date)-24*60*30,
                              max(daily_activity$local_date)+24*60*30,
                              by="days")),
             size=0.5,
             colour="darkgrey")+
  geom_vline(xintercept=max(daily_activity$local_date)+24*60*30,
             size=0.5,
             colour="black")+
  scale_x_continuous(breaks = c(seq(
                              local_date_to_posix('2022-06-11'),
                              local_date_to_posix('2022-07-11'),
                              by="days")),
                     expand=c(0,0))+
    scale_y_continuous(expand=c(0,0),
                       breaks = seq(0,45,by=10))+
    theme(axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank(),
          text=element_text(family=font.family),
          legend.title = element_text(size=legend.heading.size),
          legend.text = element_text(size=label.size),
          axis.title = element_text(size=axis.heading.size),
          axis.text = element_text(size=label.size)) +
  ylab("\n number of fertile females")

season_act_plot <- boxplot_active_date +
  hor_barplot_ff +
  plot_layout(
  nrow=1,
  widths = c(1, 0.5))

season_act_plot

Model: Male % active and number of fertile females

ff_m <- glmmTMB(prop_activity ~ 
                  scale(n_FertileFemales) + 
                  (scale(n_FertileFemales)|bird_ID),
                  data = daily_activity,
                  family = beta_family(link="logit"))
summary(ff_m)
##  Family: beta  ( logit )
## Formula:          
## prop_activity ~ scale(n_FertileFemales) + (scale(n_FertileFemales) |  
##     bird_ID)
## Data: daily_activity
## 
##      AIC      BIC   logLik deviance df.resid 
##     -590     -562      301     -602      796 
## 
## Random effects:
## 
## Conditional model:
##  Groups  Name                    Variance Std.Dev. Corr 
##  bird_ID (Intercept)             0.186    0.432         
##          scale(n_FertileFemales) 0.176    0.419    0.77 
## Number of obs: 802, groups:  bird_ID, 81
## 
## Dispersion parameter for beta family (): 7.61 
## 
## Conditional model:
##                         Estimate Std. Error z value          Pr(>|z|)    
## (Intercept)               0.2660     0.0562    4.73 0.000002228925850 ***
## scale(n_FertileFemales)   0.4326     0.0564    7.67 0.000000000000018 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(ff_m, level=0.95)
##                                                 2.5 % 97.5 % Estimate
## (Intercept)                                     0.156  0.376    0.266
## scale(n_FertileFemales)                         0.322  0.543    0.433
## Std.Dev.(Intercept)|bird_ID                     0.350  0.533    0.432
## Std.Dev.scale(n_FertileFemales)|bird_ID         0.335  0.524    0.419
## Cor.scale(n_FertileFemales).(Intercept)|bird_ID 0.471  0.886    0.774
simulationOutput <- simulateResiduals(ff_m)
testResiduals(simulationOutput)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.06, p-value = 0.01
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9, p-value = 0.04
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 10, observations = 802, p-value = 0.2
## alternative hypothesis: true probability of success is not equal to 0.00797
## 95 percent confidence interval:
##  0.0060 0.0228
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                 0.0125
ff_flex <- as_flextable(ff_m)
if(save_tables==TRUE){
save_as_docx(ff_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/fertile_females_model.docx"))
}

Male activity (% time) was higher when there were more fertile females at the site.

We can also look at a variation of this model to better understand within- and among-individual variation.

Model: Within- and among-individual variation in male activity

mean_fertile_females <- daily_activity %>%
  group_by(colour_combo, bird_ID) %>%
  summarise(ind_mean_ff = mean(n_FertileFemales))

daily_activity <- left_join(daily_activity, mean_fertile_females,
            join_by('colour_combo','bird_ID')) %>%
  mutate(ind_dev_ff=n_FertileFemales-ind_mean_ff)

In our study, there are two proximate mechanisms affecting the relationship between number of fertile females and male activity: individual changes in activity (within-individual variation) and differences in the presence of males at the site (between-individual variation). One method for distinguishing between these two mechanisms is within-subject centering (van de Pol and Wright 2009; Dingemanse and Dochtermann 2013).

From Dingemanse and Dochtermann 2013: “[The] conflation of between- and within-individual effects can be separated by calculating the mean covariate value for the individual, and the deviation of the covariate from the mean for an individual for each measurement. Both of these predictor variables are then included in modelling the phenotypic response”.

The effect of the latter (the “deviation” term) tells us the effect of the number of fertile females on male activity. The former (the “individual mean” term) tells us the effect of individuals being measured when there are different numbers of fertile females.

We can again include a random slope, so that we can also look at 1) among-individual variation in the response to changes in the number of fertile females, and 2) how activity ‘elevation’ relates to activity ‘slope’ (covariance between random intercept and slope).

ff_m2 <- glmmTMB(prop_activity ~ ind_mean_ff + ind_dev_ff +
                  (ind_dev_ff|bird_ID),
                  data = daily_activity,
                  family = beta_family(link="logit"))
summary(ff_m2)
##  Family: beta  ( logit )
## Formula:          
## prop_activity ~ ind_mean_ff + ind_dev_ff + (ind_dev_ff | bird_ID)
## Data: daily_activity
## 
##      AIC      BIC   logLik deviance df.resid 
##     -591     -558      302     -605      795 
## 
## Random effects:
## 
## Conditional model:
##  Groups  Name        Variance Std.Dev. Corr 
##  bird_ID (Intercept) 0.246    0.4956        
##          ind_dev_ff  0.001    0.0317   0.79 
## Number of obs: 802, groups:  bird_ID, 81
## 
## Dispersion parameter for beta family (): 7.71 
## 
## Conditional model:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.21055    0.20567    1.02     0.31    
## ind_mean_ff  0.00417    0.00786    0.53     0.60    
## ind_dev_ff   0.03480    0.00446    7.80    6e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(ff_m2, level=0.95)
##                                      2.5 % 97.5 % Estimate
## (Intercept)                        -0.1926 0.6137  0.21055
## ind_mean_ff                        -0.0112 0.0196  0.00417
## ind_dev_ff                          0.0261 0.0435  0.03480
## Std.Dev.(Intercept)|bird_ID         0.4054 0.6060  0.49564
## Std.Dev.ind_dev_ff|bird_ID          0.0249 0.0403  0.03168
## Cor.ind_dev_ff.(Intercept)|bird_ID  0.5369 0.8885  0.78956
simulationOutput <- simulateResiduals(ff_m)
testResiduals(simulationOutput)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.06, p-value = 0.01
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.9, p-value = 0.04
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 10, observations = 802, p-value = 0.2
## alternative hypothesis: true probability of success is not equal to 0.00797
## 95 percent confidence interval:
##  0.0060 0.0228
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                 0.0125

Based on this model, we can make the following conclusions:

  • While both within- and between-individual variation contribute to the relationship between mating opportunities and male activity, within-individual variation has a stronger effect

  • There does not seem to be that much among-individual variation in activity slopes (i.e. most males are increasing their activity by a similar amount as number of fertile females increases?)

  • There is a fairly strong positive correlation between individual ‘elevation’ and ‘slope’. So, males who were more active at their ‘centred’ number of fertile females are changing their activity more than other males.

ff2_flex <- as_flextable(ff_m2)
if(save_tables==TRUE){
save_as_docx(ff2_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/fertile_females_variation_model.docx"))
}

Fig: Male % active and number of fertile females

m_effects = effects::allEffects(ff_m, xlevels = 160)
m_ff_effects = m_effects$`scale(n_FertileFemales)` %>% data.frame
activity_ff_plot <- ggplot(daily_activity, aes(x=n_FertileFemales,
                                               y=perc_active))+
  theme_classic()+
  #geom_boxplot(colour="gray30", outlier.shape = NA)+
  geom_point(size=small.point.size,
             shape=default_shape,
             alpha=0.15,
             colour=male_colour,
             fill=male_colour)+
  geom_line(data = m_ff_effects,
            aes(y = fit*100, x = n_FertileFemales),
            colour=col_scale[1]) +
  geom_ribbon(data = m_ff_effects,
              aes(y = fit*100, x = n_FertileFemales, ymin = lower*100, ymax = upper*100),
              alpha = 0.2, fill=col_scale[1])+
  coord_cartesian(ylim = c(0,102)) +
  scale_y_continuous(breaks=seq(0,100,by=20),
                     expand = c(0,0)) +
  xlab("\nnumber of fertile females\n")+
  ylab("daily activity (% time)\n")+
  theme(axis.ticks.length.y = unit(0.07, "cm"),
        text=element_text(family=font.family),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size))

activity_ff_plot

if(save_plots==TRUE){
  ggsave(plot=activity_ff_plot,
       filename='figure-activity-fertility-means-plot.jpg',
       path=output_path,
       width=16,
       height=12.5,
       units="cm",
       dpi=300)
}

Combined figure: Date and male activity

layout <- "
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
"

fig1 <- activity_fert_boxplot + ylab("\nmale activity (% time)\n") + xlab("\nbreeding stage\n\n") +
  plot_spacer() +
  #plot_spacer() +
  season_act_plot +
  plot_layout(design=layout) +
  plot_layout(nrow=2)+
  plot_annotation(tag_levels = 'a',
                  tag_prefix = '(',
                  tag_suffix = ')') &
  theme(plot.tag = element_text(size=label.size, family = font.family, face="italic"))


fig1

if(save_plots==TRUE){
ggsave(plot=fig1,
       filename='figure1-activity-stage_date.jpg',
       path=output_path,
       width=16,
       height=20,
       units="cm",
       dpi=300)
}

Q: Do males increase the intensity of their activity when there are more mating opportunities?

daily_activity %>%
  ungroup() %>%
  filter(n_FertileFemales==0 | n_FertileFemales==40) %>%
  select(colour_combo, mean_ODBA, n_FertileFemales) %>%
  group_by(colour_combo, n_FertileFemales) %>%
  summarise(mean_ODBA=mean(mean_ODBA)) %>%
  pivot_wider(values_from = 'mean_ODBA', names_from = 'n_FertileFemales', names_prefix="ff_") %>%
  filter(!is.na(ff_40) & !is.na(ff_0)) %>%
  mutate(diff=ff_40-ff_0) %>%
  flextable()

colour_combo

ff_40

ff_0

diff

DB,O,DG

0.422

0.201

0.22086

DB,PI,PI

0.193

0.189

0.00373

DB,R,O

0.566

0.138

0.42849

DG,R,DG

0.211

0.188

0.02260

O,DB,DG

0.530

0.211

0.31961

O,DG,DG

0.505

0.180

0.32515

O,PI,Y

0.164

0.182

-0.01871

PI,DG,DG

0.596

0.205

0.39034

PI,PI,PI

0.552

0.160

0.39183

PI,Y,R

0.193

0.175

0.01812

R,DB,R

0.390

0.240

0.15035

R,Y,R

0.315

0.117

0.19782

Y,DB,Y

0.449

0.193

0.25670

Y,DG,PI

0.248

0.177

0.07112

Y,O,R

0.501

0.175

0.32574

Y,Y,DG

0.498

0.204

0.29401

The table above shows the mean ODBA of males (each with a unique colour_combo) when female fertility was at its peak (40 fertile females, ‘ff_40’) compared with when females were no longer fertile and were incubating (0 fertile females, ‘ff_0’). Almost all males showed more intense activity during the peak fertile period.

The table below summarises mean ODBA levels when there were no fertile females.

daily_activity %>%
  ungroup() %>%
  filter(n_FertileFemales==0) %>%
  summarise(sd_ODBA=sd(mean_ODBA),
            mean_ODBA=mean(mean_ODBA)) %>%
  flextable()

sd_ODBA

mean_ODBA

0.0343

0.186

The table below summarises mean ODBA levels when there were the maximum number of fertile females.

daily_activity %>%
  ungroup() %>%
  filter(n_FertileFemales==40) %>%
  summarise(sd_ODBA=sd(mean_ODBA),
            mean_ODBA=mean(mean_ODBA)) %>%
  flextable()

sd_ODBA

mean_ODBA

0.168

0.388

Prediction: Males increase the intensity of their activity (mean ODBA) when there are fertile females at the site.

Fig: Breeding stage and male ODBA

odba_fert_boxplot <- ggplot() +
  theme_classic()+
  geom_boxplot(data = ind_activity_plotdat %>% 
                 filter(plot_factor == "fertile box"),
               aes(x = plot_factor,
                   y = mean_ODBA),
               outlier.shape=default_shape,
               fill="white",
               width=0.7) +
  geom_point(data = ind_activity_plotdat %>% 
                 filter(plot_factor == "fertile point" |
                          plot_factor == "post-fertile point"),
               aes(x = plot_factor,
                   y = mean_ODBA,
                   fill=success, colour=success, alpha=success),
             size=small.point.size,
             position=position_jitter(w=0.15, h=0, seed = 1)) +
  geom_boxplot(data = ind_activity_plotdat %>% 
                 filter(plot_factor == "post-fertile box"),
               aes(x = plot_factor,
                   y = mean_ODBA),
               outlier.shape=default_shape,
               fill="white",
               width=0.7)+
  geom_line(data = ind_activity_plotdat %>% 
                 filter(plot_factor == "fertile point" |
                          plot_factor == "post-fertile point"),
            aes(x = plot_factor,
                y = mean_ODBA,
                colour=success, group=bird_ID, alpha=success), 
            size=small.point.size/3,
            position=position_jitter(w=0.15, h=0, seed = 1)) +
  coord_cartesian(ylim=c(0,0.6))+
  scale_y_continuous(breaks=c(seq(0,0.6,by=0.1)),
                     expand=c(0,0))+
  scale_x_discrete(breaks = c("fertile box", "post-fertile box"),
                   labels = c("fertile","post-fertile"))+
  scale_colour_manual(values=c(fail_colour,success_colour), name="Sired offspring?")+
  scale_fill_manual(values=c(fail_colour,success_colour), name="Sired offspring?")+
  scale_alpha_manual(values=c(0.5,0.8), name="Sired offspring?")+
  scale_shape_manual(values=c(default_shape), guide="none") +
  xlab("\n breeding stage")+
  ylab("mean male ODBA (g)\n")+
  theme(text=element_text(family=font.family),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size),
        legend.title = element_text(size=legend.heading.size),
        legend.text = element_text(size=label.size),
        legend.position="right",
        legend.justification="top",
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "grey"))

odba_fert_boxplot

ind_activity %>%
  group_by(breeding_stage, success) %>%
  summarise(mean_mean_ODBA = mean(mean_ODBA),
            sd_ODBA = sd(mean_ODBA),
            n=n(),
            se_ODBA = sd_ODBA/sqrt(n)) %>%
  flextable()

breeding_stage

success

mean_mean_ODBA

sd_ODBA

n

se_ODBA

fertile

no

0.284

0.0964

65

0.01196

fertile

yes

0.354

0.0725

16

0.01811

post-fertile

no

0.201

0.0347

41

0.00542

post-fertile

yes

0.222

0.0484

15

0.01248

In the plot above, there are 81 individual males in total, PI,PI,R, O,DB,DG, PI,Y,R, Y,DB,Y, PI,PI,O, O,DB,O, DG,O,DG, DG,DG,R, DG,R,DB, R,PI,PI, Y,Y,PI, Y,Y,O, DG,O,O, Y,DG,PI, R,O,O, DB,DG,Y, O,R,PI, R,DG,Y, DG,PI,O, PI,R,PI, DB,O,DG, PI,DG,DG, Y,Y,R, Y,O,Y, Y,Y,DG, DB,Y,PI, R,DG,R, R,O,R, PI,O,R, R,R,DG, Y,O,R, DB,Y,R, R,DG,PI, O,DG,DG, DB,R,O, DG,R,DG, Y,Y,Y, DG,R,O, O,R,O, O,PI,Y, PI,R,DB, O,DB,R, DG,DB,Y, DB,DG,DG, Y,DG,DG, DB,R,R, R,DB,R, DB,PI,PI, DG,PI,Y, R,Y,Y, PI,PI,PI, DG,R,R, DB,R,Y, R,Y,R, O,Y,O, PI,DB,DB of which have post-fertile data. Males with “uncertain” success are treated as “no” success.

Each datapoint represents the mean ODBA of an individual male.

Model: Breeding stage and male ODBA

odba_stage_m <- glmmTMB(mean_ODBA ~ breeding_stage + (breeding_stage|bird_ID),
                   data=daily_activity,
                   #data=daily_activity %>% filter(likely_undetected_success=="no"),
                   family=beta_family(link="logit"))
summary(odba_stage_m)
##  Family: beta  ( logit )
## Formula:          mean_ODBA ~ breeding_stage + (breeding_stage | bird_ID)
## Data: daily_activity
## 
##      AIC      BIC   logLik deviance df.resid 
##    -1690    -1662      851    -1702      796 
## 
## Random effects:
## 
## Conditional model:
##  Groups  Name                       Variance Std.Dev. Corr  
##  bird_ID (Intercept)                0.185    0.430          
##          breeding_stagepost-fertile 0.111    0.333    -0.99 
## Number of obs: 802, groups:  bird_ID, 81
## 
## Dispersion parameter for beta family (): 32.4 
## 
## Conditional model:
##                            Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                 -0.8712     0.0511  -17.03   <2e-16 ***
## breeding_stagepost-fertile  -0.4526     0.0493   -9.18   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(odba_stage_m, level=0.95)
##                                                     2.5 % 97.5 % Estimate
## (Intercept)                                        -0.971 -0.771   -0.871
## breeding_stagepost-fertile                         -0.549 -0.356   -0.453
## Std.Dev.(Intercept)|bird_ID                         0.361  0.512    0.430
## Std.Dev.breeding_stagepost-fertile|bird_ID          0.258  0.431    0.333
## Cor.breeding_stagepost-fertile.(Intercept)|bird_ID -0.999  0.991   -0.986
simulationOutput <- simulateResiduals(odba_stage_m)
testResiduals(simulationOutput)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.05, p-value = 0.02
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1, p-value = 0.1
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 6, observations = 802, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.00797
## 95 percent confidence interval:
##  0.00275 0.01621
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                0.00748
odba_stage_m_flex <- as_flextable(odba_stage_m)
odba_stage_m_flex

group

Estimate

Standard Error

statistic

p-value

Fixed effects

(Intercept)

-0.871

0.051

-17.034

0.0000

***

breeding_stagepost-fertile

-0.453

0.049

-9.175

0.0000

***

Random effects

bird_ID

sd__(Intercept)

0.430

bird_ID

sd__breeding_stagepost-fertile

0.333

bird_ID

cor__(Intercept).breeding_stagepost-fertile

-0.986

Signif. codes: 0 <= '***' < 0.001 < '**' < 0.01 < '*' < 0.05

square root of the estimated residual variance: 32.4

data's log-likelihood under the model: 850.8

Akaike Information Criterion: -1,689.7

Bayesian Information Criterion: -1,661.6

if(save_tables==TRUE){
save_as_docx(odba_stage_m_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/breeding_stage_odba_model.docx"))
}

When ‘potentially’ successful males are excluded the conclusions remain the same.

as.data.frame(emmeans(odba_stage_m, ~ breeding_stage))
##  breeding_stage emmean     SE  df asymp.LCL asymp.UCL
##  fertile        -0.871 0.0511 Inf    -0.971    -0.771
##  post-fertile   -1.324 0.0307 Inf    -1.384    -1.264
## 
## Results are given on the logit (not the response) scale. 
## Confidence level used: 0.95
odba_stage_emm_flex <- as_flextable(as.data.frame(emmeans(odba_stage_m, ~ breeding_stage)))

if(save_tables==TRUE){
save_as_docx(odba_stage_emm_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/breeding_stage_emm.docx"))
}

Male tended to show more intense activity during the fertile stage, compared with the post-fertile stage. The inclusion of males that were more likely to have sired additional offspring (predated before they could be sampled) had no effect on these results.

Prediction 2 Intensity of male activity (mean ODBA) increases with the number of fertile females at the site.

We also examine this trend in more detail, by looking at activity and the number of fertile females each day across the season.

Fig: Date and male ODBA

boxplot_odba_date <- ggplot(daily_activity, aes(x=local_date,
                                                 y=mean_ODBA,
                                                 group=local_date))+
  geom_boxplot(colour="gray30", outlier.shape = NA)+
  geom_point(size=small.point.size,
             aes(colour=success, alpha=success, fill=success),
             shape=default_shape)+
  geom_vline(xintercept=c(seq(min(daily_activity$local_date)-24*60*30,
                              max(daily_activity$local_date)+24*60*30,
                              by="days")),
             size=0.5,
             colour="darkgrey")+
  geom_vline(xintercept=max(daily_activity$local_date)+24*60*30,
             size=0.5,
             colour="black")+
  coord_flip(ylim=c(0,0.9),
             xlim = c(max(daily_activity$local_date)+24*60*30,
                      min(daily_activity$local_date)-24*60*30))+
  scale_x_continuous(breaks = c(seq(
                              local_date_to_posix('2022-06-11'),
                              local_date_to_posix('2022-07-11'),
                              by="days")),
                     labels=date_labels,
                     expand=c(0,0))+
  scale_y_continuous(breaks=seq(0,0.9,by=0.2),
                     expand=c(0,0)) +
  scale_colour_manual(values = c(fail_colour, success_colour), guide="none")+
  scale_fill_manual(values = c(fail_colour, success_colour), guide="none")+
  scale_alpha_manual(values = c(0.5, 0.8), guide="none")+
  xlab("date\n")+
  ylab("\nmean male ODBA (g)")+
  geom_text(data = date_activity_n,
            family=font.family,
            size=(label.size-2)/.pt,
            colour="gray50",
            aes(local_date, 107, label = n_ind),
            vjust = 0.3,
            hjust = 0)+
  theme(axis.ticks.length.y = unit(0.07, "cm"),
        text=element_text(family=font.family),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size),
        axis.line.y = element_line(size=0.5, color="black"),
        panel.background=element_blank())

boxplot_odba_date

season_odba_plot <- boxplot_odba_date +
  hor_barplot_ff +
  plot_layout(
  nrow=1,
  widths = c(1, 0.5))

season_odba_plot

Model: Fertile females and ODBA

ff_m_odba <- glmmTMB(mean_ODBA ~ scale(n_FertileFemales) +
                     (scale(n_FertileFemales)|bird_ID),
                  data = daily_activity,
                  family = beta_family(link="logit"))
summary(ff_m_odba)
##  Family: beta  ( logit )
## Formula:          
## mean_ODBA ~ scale(n_FertileFemales) + (scale(n_FertileFemales) |      bird_ID)
## Data: daily_activity
## 
##      AIC      BIC   logLik deviance df.resid 
##    -1764    -1736      888    -1776      796 
## 
## Random effects:
## 
## Conditional model:
##  Groups  Name                    Variance Std.Dev. Corr 
##  bird_ID (Intercept)             0.1269   0.356         
##          scale(n_FertileFemales) 0.0412   0.203    0.63 
## Number of obs: 802, groups:  bird_ID, 81
## 
## Dispersion parameter for beta family (): 39.1 
## 
## Conditional model:
##                         Estimate Std. Error z value Pr(>|z|)    
## (Intercept)              -1.0219     0.0429  -23.80   <2e-16 ***
## scale(n_FertileFemales)   0.2398     0.0287    8.36   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(ff_m_odba, level=0.95)
##                                                  2.5 % 97.5 % Estimate
## (Intercept)                                     -1.106 -0.938   -1.022
## scale(n_FertileFemales)                          0.184  0.296    0.240
## Std.Dev.(Intercept)|bird_ID                      0.298  0.426    0.356
## Std.Dev.scale(n_FertileFemales)|bird_ID          0.159  0.259    0.203
## Cor.scale(n_FertileFemales).(Intercept)|bird_ID  0.318  0.788    0.628
simulationOutput <- simulateResiduals(ff_m_odba)
testResiduals(simulationOutput)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.06, p-value = 0.008
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1, p-value = 0.2
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 10, observations = 802, p-value = 0.2
## alternative hypothesis: true probability of success is not equal to 0.00797
## 95 percent confidence interval:
##  0.0060 0.0228
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                 0.0125
ff_odba_flex <- as_flextable(ff_m_odba)
if(save_tables==TRUE){
save_as_docx(ff_odba_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/fertile_females_odba_model.docx"))
}

Model: Within- and among-individual variation

ff_m_odba2 <- glmmTMB(mean_ODBA ~ ind_mean_ff + ind_dev_ff +
                  (ind_dev_ff|bird_ID),
                  data = daily_activity,
                  family = beta_family(link="logit"))
summary(ff_m_odba2)
##  Family: beta  ( logit )
## Formula:          mean_ODBA ~ ind_mean_ff + ind_dev_ff + (ind_dev_ff | bird_ID)
## Data: daily_activity
## 
##      AIC      BIC   logLik deviance df.resid 
##    -1768    -1736      891    -1782      795 
## 
## Random effects:
## 
## Conditional model:
##  Groups  Name        Variance Std.Dev. Corr 
##  bird_ID (Intercept) 0.143905 0.379         
##          ind_dev_ff  0.000226 0.015    0.68 
## Number of obs: 802, groups:  bird_ID, 81
## 
## Dispersion parameter for beta family (): 39.5 
## 
## Conditional model:
##             Estimate Std. Error z value       Pr(>|z|)    
## (Intercept) -1.03498    0.15620   -6.63 0.000000000034 ***
## ind_mean_ff  0.00197    0.00592    0.33           0.74    
## ind_dev_ff   0.01817    0.00218    8.34        < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(ff_m_odba2, level=0.95)
##                                       2.5 %  97.5 % Estimate
## (Intercept)                        -1.34113 -0.7288 -1.03498
## ind_mean_ff                        -0.00964  0.0136  0.00197
## ind_dev_ff                          0.01391  0.0224  0.01817
## Std.Dev.(Intercept)|bird_ID         0.31889  0.4513  0.37935
## Std.Dev.ind_dev_ff|bird_ID          0.01172  0.0193  0.01504
## Cor.ind_dev_ff.(Intercept)|bird_ID  0.41008  0.8132  0.67846
simulationOutput <- simulateResiduals(ff_m_odba2)
testResiduals(simulationOutput)

## $uniformity
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.05, p-value = 0.04
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1, p-value = 0.4
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa outlier test based on exact binomial test with approximate
##  expectations
## 
## data:  simulationOutput
## outliers at both margin(s) = 10, observations = 802, p-value = 0.2
## alternative hypothesis: true probability of success is not equal to 0.00797
## 95 percent confidence interval:
##  0.0060 0.0228
## sample estimates:
## frequency of outliers (expected: 0.00796812749003984 ) 
##                                                 0.0125

Based on this model, we can make the following conclusions:

  • The relationship between mating opportunities and male activity is explained by within-individual variation

  • There does not seem to be that much among-individual variation in activity slopes (i.e. most males are increasing their activity by a similar amount as number of fertile females increases?)

  • There is a fairly strong negative correlation between individual ‘elevation’ and ‘slope’, like what we saw with percent time active

ff_odba2_flex <- as_flextable(ff_m_odba2)
if(save_tables==TRUE){
save_as_docx(ff_odba2_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/fertile_females_odba_variation_model.docx"))
}

Fig: Male ODBA and number of fertile females

m_effects = effects::allEffects(ff_m_odba, xlevels = 160)
m_ff_odba_effects = m_effects$`scale(n_FertileFemales)` %>% data.frame

odba_ff_plot <- ggplot(daily_activity, aes(x=n_FertileFemales,
                                               y=mean_ODBA))+
  theme_classic()+
  #geom_boxplot(colour="gray30", outlier.shape = NA)+
  geom_point(size=small.point.size,
             shape=default_shape,
             alpha=0.15,
             colour=male_colour,
             fill=male_colour)+
  geom_line(data = m_ff_odba_effects,
            aes(y = fit, x = n_FertileFemales),
            colour=col_scale[1]) +
  geom_ribbon(data = m_ff_odba_effects,
              aes(y = fit, x = n_FertileFemales, ymin = lower, ymax = upper),
              alpha = 0.2, fill=col_scale[1])+
  coord_cartesian(ylim = c(0,0.9)) +
  scale_y_continuous(breaks=seq(0,0.9,by=0.2),
                     expand = c(0,0)) +
  xlab("\nnumber of fertile females\n")+
  ylab("mean daily ODBA (g)\n")+
  theme(axis.ticks.length.y = unit(0.07, "cm"),
        text=element_text(family=font.family),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size))

odba_ff_plot

daily_activity %>%
  ungroup() %>%
  filter(n_FertileFemales == 0 | n_FertileFemales==40) %>%
  group_by(n_FertileFemales) %>%
  summarise(mean_mean_ODBA = mean(mean_ODBA),
            sd_odba_active = sd(mean_ODBA),
            n=n(),
            se_odba_active=sd_odba_active/sqrt(n))
## # A tibble: 2 × 5
##   n_FertileFemales mean_mean_ODBA sd_odba_active     n se_odba_active
##              <int>          <dbl>          <dbl> <int>          <dbl>
## 1                0          0.186         0.0343    73        0.00402
## 2               40          0.388         0.168    122        0.0152

Combined figure: Date and male ODBA

layout <- "
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
AAAAAAAAAAAAAAAAAAAABB
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
DDDDDDDDDDDDDDDDDDDDDD
"

fig_s1 <- odba_fert_boxplot + ylab("\nmean male ODBA (g)\n") + xlab("\nbreeding stage\n\n") +
  plot_spacer() +
  #plot_spacer() +
  season_odba_plot +
  plot_layout(design=layout) +
  plot_layout(nrow=2)+
  plot_annotation(tag_levels = 'a',
                  tag_prefix = '(',
                  tag_suffix = ')') &
  theme(plot.tag = element_text(size=label.size, family = font.family, face="italic"))


fig_s1

if(save_plots==TRUE){
ggsave(plot=fig_s1,
       filename='figure-odba-stage_date.jpg',
       path=output_path,
       width=16,
       height=20,
       units="cm",
       dpi=300)
}

Combined figure: Fertile females and male activity

fig2 <- activity_ff_plot +
  plot_spacer() +
  odba_ff_plot +
  plot_layout(nrow = 1,
              widths=c(1,0.1,1),
              heights=c(1))+
  plot_annotation(tag_levels = 'a',
                  tag_prefix = '(',
                  tag_suffix = ')') &
  theme(plot.tag = element_text(size=label.size, family = font.family, face="italic"),
        plot.tag.position = c(.27, .96))

fig2

if(save_plots==TRUE){
ggsave(plot=fig2,
       filename='figure2-activity-fertilefemales.jpg',
       path=output_path,
       width=16,
       height=7,
       units="cm",
       dpi=300)
}
rm(stage_m, stage_m_flex,
   m_ff_odba_effects, m_ff_effects, m_effects,
   ff_flex, ff_m, ff_m_odba, ff_m2,
   ff_odba_flex, ff2_flex)

Q: How does male behaviour relate to mating and reproductive success?

Prediction Increased time spent active (%) is associated with increased likelihood of mating success Prediction Males who spend more days at the site during the fertile period have a higher probability of mating success Prediction Males with higher mean ODBA have a higher probability of mating success.

In the plots and analyses below, we only include data from the fertile stage, and only for males with at least 2 days of data during the fertile stage.

However, it’s a problem to exclude males with limited data for the tenure analysis, because then we’re excluding all the males who stayed for less than 2. So we analyse the effects of tenure on mating success in a separate model, including the full dataset.

ind_activity <- ind_activity %>%
  mutate(prop_activity=total_perc_active/100,
         recording_hours=n_windows*window_length/60,
         arrival_day=as.numeric(difftime(arrival_date, 
                                         local_date_to_posix("2022-06-01"),
                                         units="days")),
         success_num = case_when(success=="yes" ~ 1,
                                 success=="no" ~ 0))

ind_fertile_activity <- ind_activity %>%
  filter(breeding_stage=="fertile" & n_days>=min_acc_days)

ind_tenure <- left_join(male_tenure, father_summary, 
                        join_by('colour_combo'=='LL.father',
                                'bird_ID'=='ID.father')) %>%
    mutate(arrival_day=as.numeric(difftime(arrival_date, 
                                         local_date_to_posix("2022-06-01"),
                                         units="days")),
           likely_extra_nests=case_when(potential_n_extra_nests>0 ~ "yes",
                                      TRUE ~ "less likely"),) %>%
    replace(is.na(.), 0)
ind_fertile_activity %>%
  ungroup() %>%
  summarise(min_total_perc_active=min(total_perc_active),
            max_total_perc_active=max(total_perc_active),
            mean_total_perc_active=mean(total_perc_active),
            sd_total_perc_active=sd(total_perc_active),
            n=n(),
            se_total_perc_active=sd_total_perc_active/sqrt(n)) %>%
  flextable()

min_total_perc_active

max_total_perc_active

mean_total_perc_active

sd_total_perc_active

n

se_total_perc_active

15.4

84.5

59.9

14.8

81

1.64

The values below show the correlation between variables considered for models.

cor(ind_fertile_activity$tenure_days, ind_fertile_activity$recording_hours)
## [1] 0.749
cor(ind_fertile_activity$tenure_days, ind_fertile_activity$prop_activity)
## [1] 0.19
cor(ind_fertile_activity$prop_activity, ind_fertile_activity$arrival_day)
## [1] -0.228
cor(ind_fertile_activity$tenure_days, ind_fertile_activity$arrival_day)
## [1] -0.201

Tenure is highly correlated with number of hours of data recorded, which makes sense. So we do not include both these variables in the model. Tenure is the more interesting and biologically meaningful variable, so we focus on that.

Also note that mean ODBA and proportion of time spent active are correlated.

cor.test(ind_fertile_activity$mean_ODBA, ind_fertile_activity$prop_activity)
## 
##  Pearson's product-moment correlation
## 
## data:  ind_fertile_activity$mean_ODBA and ind_fertile_activity$prop_activity
## t = 19, df = 79, p-value <2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.853 0.937
## sample estimates:
##   cor 
## 0.903

Mating success (number of females with which a male sired offspring)

Model: % active and mating success

There are often concerns about models with Poisson distributions. We confirmed that all model assumptions were met and that there were no issues with overdispersion. Running the model with a Gaussian distribution violates model assumptions but still produces the same conclusions.

success_m1 <- glmmTMB(n_nests ~ scale(prop_activity) + scale(log(tenure_days)) + scale(arrival_day),
                   data=ind_fertile_activity,
                   # data=ind_fertile_activity %>% filter(likely_undetected_success=="no"),
                   family=poisson)
summary(success_m1)
##  Family: poisson  ( log )
## Formula:          
## n_nests ~ scale(prop_activity) + scale(log(tenure_days)) + scale(arrival_day)
## Data: ind_fertile_activity
## 
##      AIC      BIC   logLik deviance df.resid 
##     97.8    107.4    -44.9     89.8       77 
## 
## 
## Conditional model:
##                         Estimate Std. Error z value     Pr(>|z|)    
## (Intercept)               -1.694      0.293   -5.78 0.0000000073 ***
## scale(prop_activity)       0.436      0.264    1.66        0.098 .  
## scale(log(tenure_days))    0.554      0.332    1.67        0.095 .  
## scale(arrival_day)         0.109      0.253    0.43        0.666    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(success_m1, level=0.95)
##                           2.5 % 97.5 % Estimate
## (Intercept)             -2.2682 -1.120   -1.694
## scale(prop_activity)    -0.0802  0.953    0.436
## scale(log(tenure_days)) -0.0956  1.204    0.554
## scale(arrival_day)      -0.3874  0.606    0.109
simulationOutput <- simulateResiduals(fittedModel = success_m1, plot=F)
testResiduals(simulationOutput)

## $uniformity
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.07, p-value = 0.8
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1, p-value = 0.7
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 81, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0000 0.0123
## sample estimates:
## outlier frequency (expected: 0.00197530864197531 ) 
##                                                  0
success_m1_flex <- as_flextable(success_m1)
if(save_tables==TRUE){
save_as_docx(success_m1_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/activity_mating_model.docx"))
}

The effect of activity is similar when the males with territories around predated (unsampled) nests are excluded, although effect of tenure becomes larger.

m1_effects = effects::allEffects(success_m1, xlevels = 100)
m1_active_eff = m1_effects$`scale(prop_activity)` %>% data.frame
m1_tenure_eff = m1_effects$`scale(log(tenure_days))` %>% data.frame

m1_effects2 = effects::allEffects(success_m1, xlevels = 100, confidence.level=0.95)
m1_active_eff2 = m1_effects2$`scale(prop_activity)` %>% data.frame
m1_effects = as.data.frame(confint(success_m1, level=0.95))

write.csv(m1_effects,
          file=glue("{output_supp_effects}/activity-mating-effects_{window_length}min_{data_type}-{threshold_type}thresh.csv"))

Fig: % active and mating success

success_active_plot <- ggplot(ind_fertile_activity)+
  geom_line(data = m1_active_eff,
            aes(y = fit, x = prop_activity*100),
            colour=col_scale[2]) +
  geom_ribbon(data = m1_active_eff,
              aes(y = fit, x =prop_activity*100, ymin = lower, ymax = upper),
              alpha = 0.2, fill=col_scale[2])+
  geom_jitter(aes(x=prop_activity*100, y=n_nests, shape=likely_extra_nests),
             fill=male_colour,
             colour=male_colour,
             size=point_size,
             alpha=point_alpha,
             width=0.005, height=0.04) +
  theme_classic()+
  scale_shape_manual(values=c(default_shape, default_shape), guide="none")+
  scale_x_continuous(expand = c(0, 0), breaks=c(seq(0,100,by=20))) +
  scale_y_continuous(breaks=c(seq(0,2,by=1))) +
  coord_cartesian(xlim=c(10,90), ylim=c(0,2.1))+
  xlab("\n activity (% time)") +
  ylab("mating success\n")+
  theme(legend.position="right",
        plot.margin=unit(c(0,0,0,0), 'lines'),
        text=element_text(family=font.family),
        legend.title = element_text(size=legend.heading.size),
        legend.text = element_text(size=label.size),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size))

success_active_plot

Model: Tenure and mating success

success_mt <- glmmTMB(n_nests ~ scale(log(tenure_days_exact)) + scale(arrival_day),
                   data=ind_tenure,
                   # data=ind_tenure %>% filter(potential_n_extra_nests==0),
                   family=poisson)

summary(success_mt)
##  Family: poisson  ( log )
## Formula:          n_nests ~ scale(log(tenure_days_exact)) + scale(arrival_day)
## Data: ind_tenure
## 
##      AIC      BIC   logLik deviance df.resid 
##     99.7    107.5    -46.9     93.7       96 
## 
## 
## Conditional model:
##                               Estimate Std. Error z value   Pr(>|z|)    
## (Intercept)                    -2.1857     0.3982   -5.49 0.00000004 ***
## scale(log(tenure_days_exact))   1.4240     0.5590    2.55      0.011 *  
## scale(arrival_day)              0.0242     0.2527    0.10      0.924    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(success_mt, level=0.95)
##                                2.5 % 97.5 % Estimate
## (Intercept)                   -2.966 -1.405  -2.1857
## scale(log(tenure_days_exact))  0.328  2.520   1.4240
## scale(arrival_day)            -0.471  0.519   0.0242
simulationOutput <- simulateResiduals(fittedModel = success_mt, plot=F)
plot(simulationOutput)

testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1, p-value = 0.7
## alternative hypothesis: two.sided
success_mt_flex <- as_flextable(success_mt)
if(save_tables==TRUE){
save_as_docx(success_mt_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/tenure_mating_model.docx"))
}

The effect of activity is similar when the males with territories around predated (unsampled) nests are excluded, although effect of tenure becomes larger.

mt_effects = effects::allEffects(success_mt, xlevels = 100)
mt_tenure_eff = mt_effects$`scale(log(tenure_days_exact))` %>% data.frame

Fig: Tenure and mating success

success_tenure_plot <- ggplot(ind_tenure)+
  geom_line(data = mt_tenure_eff,
            aes(y = fit, x = tenure_days_exact),
            colour=col_scale[1]) +
  geom_ribbon(data = mt_tenure_eff,
              aes(y = fit, x =tenure_days_exact, ymin = lower, ymax = upper),
              alpha = 0.2, fill=col_scale[1]) +
  geom_jitter(aes(x=tenure_days_exact, y=n_nests, shape=likely_extra_nests),
             fill=male_colour,
             colour=male_colour,
             size=point_size,
             alpha=point_alpha,
             width=0.01, height=0.05) +
  theme_classic()+
  scale_shape_manual(values=c(default_shape, default_shape), guide="none")+
  scale_x_continuous(expand = c(0, 0), breaks=c(seq(0,30,by=10))) +
  coord_cartesian(xlim=c(0,32), ylim=c(0,2.1))+
  scale_y_continuous(breaks=c(seq(0,2,by=1))) +
  xlab("\n tenure (days)") +
  ylab("mating success\n")+
  theme(legend.position="right",
        plot.margin=unit(c(0,0,0,0), 'lines'),
        text=element_text(family=font.family),
        legend.title = element_text(size=legend.heading.size),
        legend.text = element_text(size=label.size),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size))

success_tenure_plot

I double-checked these long tenures of unsuccessful males - they were mostly at the site or very close by, including during the fertile period. In other words, they had plenty of opportunity to sire offspring at the study site, but didn’t.

ind_tenure %>%
  ungroup() %>%
  summarise(mean_tenure=mean(tenure_days_exact),
            sd_tenure=sd(tenure_days_exact),
            n=n(),
            se_tenure=sd_tenure/sqrt(n))
##   mean_tenure sd_tenure  n se_tenure
## 1        9.49      6.25 99     0.628
ind_fertile_activity %>%
  filter(success=="yes") %>%
  select(tenure_days_exact, n_nests) %>%
  arrange(tenure_days_exact)
## # A tibble: 16 × 10
## # Groups:   bird_ID, colour_combo, breeding_stage, body_mass, n_eggs, n_nests,
## #   success, likely_undetected_success, potential_n_extra_nests [16]
##      bird_ID colour_combo breeding_stage body_mass n_eggs success
##        <int> <chr>        <chr>              <dbl>  <int> <chr>  
##  1 141269127 DG,DB,DG     fertile             96.6      4 yes    
##  2 141269182 O,Y,O        fertile             97.9      2 yes    
##  3 141269175 DG,R,R       fertile             98.5      3 yes    
##  4 141269043 DG,O,O       fertile            101.       4 yes    
##  5 141269152 Y,Y,Y        fertile            100.       1 yes    
##  6 141269146 R,DG,PI      fertile            102.       2 yes    
##  7 141269046 DB,DG,Y      fertile             95.3      4 yes    
##  8 141269169 DB,R,R       fertile             97.5      4 yes    
##  9 141269026 O,DB,O       fertile             94.7      4 yes    
## 10 141269125 Y,Y,R        fertile             95.8      1 yes    
## 11 141269144 DB,Y,R       fertile            102.       3 yes    
## 12 141269147 O,DG,DG      fertile            105.       5 yes    
## 13 141269148 DB,R,O       fertile            108.       3 yes    
## 14 141269044 Y,DG,PI      fertile             92.5      5 yes    
## 15 141269114 PI,R,PI      fertile            106        3 yes    
## 16 141269009 PI,PI,R      fertile            101.       6 yes    
## # ℹ 4 more variables: likely_undetected_success <chr>,
## #   potential_n_extra_nests <int>, tenure_days_exact <dbl>, n_nests <int>

Model: ODBA and mating success

success_m2 <- glmmTMB(n_nests ~ scale(mean_ODBA) + scale(log(tenure_days)) + scale(arrival_day),
                   data=ind_fertile_activity,
                   # data=ind_fertile_activity %>% filter(likely_undetected_success=="no"),
                   family=poisson)
                   #family=gaussian)
summary(success_m2)
##  Family: poisson  ( log )
## Formula:          
## n_nests ~ scale(mean_ODBA) + scale(log(tenure_days)) + scale(arrival_day)
## Data: ind_fertile_activity
## 
##      AIC      BIC   logLik deviance df.resid 
##     95.2    104.8    -43.6     87.2       77 
## 
## 
## Conditional model:
##                         Estimate Std. Error z value    Pr(>|z|)    
## (Intercept)               -1.786      0.319   -5.59 0.000000023 ***
## scale(mean_ODBA)           0.624      0.283    2.21       0.027 *  
## scale(log(tenure_days))    0.585      0.353    1.66       0.098 .  
## scale(arrival_day)         0.257      0.269    0.96       0.338    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(success_m2, level=0.95)
##                           2.5 % 97.5 % Estimate
## (Intercept)             -2.4121 -1.160   -1.786
## scale(mean_ODBA)         0.0697  1.179    0.624
## scale(log(tenure_days)) -0.1072  1.277    0.585
## scale(arrival_day)      -0.2692  0.784    0.257
simulationOutput <- simulateResiduals(fittedModel = success_m2, plot=F)
testResiduals(simulationOutput)

## $uniformity
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.07, p-value = 0.8
## alternative hypothesis: two-sided
## 
## 
## $dispersion
## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1, p-value = 0.8
## alternative hypothesis: two.sided
## 
## 
## $outliers
## 
##  DHARMa bootstrapped outlier test
## 
## data:  simulationOutput
## outliers at both margin(s) = 0, observations = 81, p-value = 1
## alternative hypothesis: two.sided
##  percent confidence interval:
##  0.0000 0.0123
## sample estimates:
## outlier frequency (expected: 0.00172839506172839 ) 
##                                                  0
plot(simulationOutput)

success_m2_flex <- as_flextable(success_m2)
if(save_tables==TRUE){
save_as_docx(success_m2_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/odba_mating_model.docx"))
}

Similar effects are observed when males near predated nests are excluded.

m2_effects = effects::allEffects(success_m2, xlevels = 100)
m2_active_eff = m2_effects$`scale(mean_ODBA)` %>% data.frame
m2_tenure_eff = m2_effects$`scale(log(tenure_days))` %>% data.frame

Fig: ODBA and mating success

success_odba_plot <- ggplot(ind_fertile_activity)+
  geom_line(data = m2_active_eff,
            aes(y = fit, x = mean_ODBA),
            colour=col_scale[1]) +
  geom_ribbon(data = m2_active_eff,
              aes(y = fit, x = mean_ODBA, ymin = lower, ymax = upper),
              alpha = 0.2, fill=col_scale[1])+
  geom_jitter(aes(x=mean_ODBA, y=n_nests, shape=likely_extra_nests),
             fill=male_colour,
             colour=male_colour,
             size=point_size,
             alpha=point_alpha,
             width=0.005, height=0.04) +
  theme_classic()+
  scale_shape_manual(values=c(default_shape, default_shape), guide="none")+
  scale_x_continuous(expand = c(0, 0), breaks=c(seq(0.1,0.5,by=0.1))) +
  scale_y_continuous(breaks=c(seq(0,2,by=1))) +
  coord_cartesian(xlim=c(0.1,0.53), ylim=c(0,2.1))+
  xlab("\n mean ODBA (g)") +
  ylab("mating success\n")+
  theme(legend.position="right",
        plot.margin=unit(c(0,0,0,0), 'lines'),
        text=element_text(family=font.family),
        legend.title = element_text(size=legend.heading.size),
        legend.text = element_text(size=label.size),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size))

success_odba_plot

Reproductive success (number of offspring sired)

Prediction 5 Increased time spent active (%) is associated with increased reproductive success Prediction 6 Males with longer site tenure have higher reproductive success Prediction 7 Males with longer site tenure have higher reproductive success

Model: % active and reproductive success

success_m3 <- glmmTMB(n_eggs ~ scale(prop_activity) + scale(log(tenure_days)) + scale(arrival_day),
                  data=ind_fertile_activity,
                  # data=ind_fertile_activity %>% filter(likely_extra_nests=="less likely"),
                  family=nbinom2)
summary(success_m3)
##  Family: nbinom2  ( log )
## Formula:          
## n_eggs ~ scale(prop_activity) + scale(log(tenure_days)) + scale(arrival_day)
## Data: ind_fertile_activity
## 
##      AIC      BIC   logLik deviance df.resid 
##    158.8    170.7    -74.4    148.8       76 
## 
## 
## Dispersion parameter for nbinom2 family (): 0.149 
## 
## Conditional model:
##                         Estimate Std. Error z value Pr(>|z|)  
## (Intercept)               -0.611      0.335   -1.82    0.068 .
## scale(prop_activity)       0.306      0.416    0.74    0.462  
## scale(log(tenure_days))    0.603      0.435    1.38    0.167  
## scale(arrival_day)        -0.026      0.301   -0.09    0.931  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
simulationOutput <- simulateResiduals(fittedModel = success_m3, plot=F)
plot(simulationOutput)

success_m3_flex <- as_flextable(success_m3)
if(save_tables==TRUE){
save_as_docx(success_m3_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/activity_reproduction_model.docx"))
}

Similar effects are observed when males near predated nests are excluded.

m3_effects = effects::allEffects(success_m3, xlevels = 100)
m3_active_eff = m3_effects$`scale(prop_activity)` %>% data.frame
m3_tenure_eff = m3_effects$`scale(log(tenure_days))` %>% data.frame
m3_effects = as.data.frame(confint(success_m3, level=0.95))

write.csv(m3_effects,
          file=glue("{output_supp_effects}/activity-repr-effects_{window_length}min_{data_type}-{threshold_type}thresh.csv"))

Fig: % active and reproductive success

egg_active_plot <- ggplot(ind_fertile_activity)+
  geom_line(data = m3_active_eff,
            aes(y = fit, x = prop_activity*100),
            colour=col_scale[2]) +
  geom_ribbon(data = m3_active_eff,
              aes(y = fit, x =prop_activity*100, ymin = lower, ymax = upper),
              alpha = 0.2, fill=col_scale[2])+
  geom_jitter(aes(x=prop_activity*100, y=n_eggs, shape=likely_extra_nests),
             fill=male_colour,
             colour=male_colour,
             size=point_size,
             alpha=point_alpha,
             width=0.005, height=0.15) +
  theme_classic()+
  scale_shape_manual(values=c(default_shape, default_shape), guide="none")+
  scale_x_continuous(expand = c(0, 0), breaks=c(seq(0,90,by=20))) +
  coord_cartesian(xlim=c(10,90), ylim=c(0,6.2))+
  xlab("\n activity (% time)") +
  ylab("reproductive success\n")+
  theme(legend.position="right",
        plot.margin=unit(c(0,0,0,0), 'lines'),
        text=element_text(family=font.family),
        legend.title = element_text(size=legend.heading.size),
        legend.text = element_text(size=label.size),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size))

egg_active_plot

Model: Tenure and reproductive success

success_mt2 <- glmmTMB(n_eggs ~ scale(log(tenure_days_exact)) + scale(arrival_day),
                   data=ind_tenure,
                   # data=ind_tenure %>% filter(potential_n_extra_nests==0),
                   family=nbinom2)
summary(success_mt2)
##  Family: nbinom2  ( log )
## Formula:          n_eggs ~ scale(log(tenure_days_exact)) + scale(arrival_day)
## Data: ind_tenure
## 
##      AIC      BIC   logLik deviance df.resid 
##    158.7    169.1    -75.3    150.7       95 
## 
## 
## Dispersion parameter for nbinom2 family (): 0.139 
## 
## Conditional model:
##                               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)                    -1.3452     0.4366   -3.08   0.0021 **
## scale(log(tenure_days_exact))   1.9216     0.7530    2.55   0.0107 * 
## scale(arrival_day)             -0.0294     0.3356   -0.09   0.9301   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(success_mt2, level=0.95)
##                                2.5 % 97.5 % Estimate
## (Intercept)                   -2.201 -0.489  -1.3452
## scale(log(tenure_days_exact))  0.446  3.397   1.9216
## scale(arrival_day)            -0.687  0.628  -0.0294
simulationOutput <- simulateResiduals(fittedModel = success_mt2, plot=F)
plot(simulationOutput)

testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.3, p-value = 0.4
## alternative hypothesis: two.sided

Same conclusions when males near predated nests are excluded.

success_mt2_flex <- as_flextable(success_mt2)
if(save_tables==TRUE){
save_as_docx(success_mt2_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/tenure_reproduction_model.docx"))
}
mt2_effects = effects::allEffects(success_mt2, xlevels = 100)
mt2_tenure_eff = mt2_effects$`scale(log(tenure_days_exact))` %>% data.frame

Fig: Tenure and reproductive success

egg_tenure_plot <- ggplot(ind_tenure) +
  geom_line(data = mt2_tenure_eff,
            aes(y = fit, x = tenure_days_exact),
            colour=col_scale[1]) +
  geom_ribbon(data = mt2_tenure_eff,
              aes(y = fit, x = tenure_days_exact, ymin = lower, ymax = upper),
              alpha = 0.2, fill=col_scale[1])+
  geom_jitter(aes(x=tenure_days_exact, y=n_eggs, shape=likely_extra_nests),
             fill=male_colour,
             colour=male_colour,
             size=point_size,
             alpha=point_alpha,
             width=0.3, height=0.15) +
  theme_classic()+
  scale_shape_manual(values=c(default_shape, default_shape), guide="none")+
  scale_x_continuous(expand = c(0, 0), breaks=c(seq(0,30,by=10))) +
  coord_cartesian(xlim=c(0,32), ylim=c(0,6.2))+
  xlab("\n tenure (days)") +
  ylab("reproductive success\n") +
  theme(legend.position="right",
        plot.margin=unit(c(0,0,0,0), 'lines'),
        text=element_text(family=font.family),
        legend.title = element_text(size=legend.heading.size),
        legend.text = element_text(size=label.size),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size))

egg_tenure_plot

Model: ODBA and reproductive success

success_m4 <- glmmTMB(n_eggs ~ scale(mean_ODBA) + scale(log(tenure_days)) + scale(arrival_day),
                  data=ind_fertile_activity,
                  # data=ind_fertile_activity %>% filter(likely_extra_nests=="less likely"),
                  family=nbinom2)
summary(success_m4)
##  Family: nbinom2  ( log )
## Formula:          
## n_eggs ~ scale(mean_ODBA) + scale(log(tenure_days)) + scale(arrival_day)
## Data: ind_fertile_activity
## 
##      AIC      BIC   logLik deviance df.resid 
##    157.2    169.2    -73.6    147.2       76 
## 
## 
## Dispersion parameter for nbinom2 family (): 0.163 
## 
## Conditional model:
##                         Estimate Std. Error z value Pr(>|z|)  
## (Intercept)               -0.692      0.331   -2.09    0.037 *
## scale(mean_ODBA)           0.629      0.425    1.48    0.139  
## scale(log(tenure_days))    0.483      0.418    1.16    0.247  
## scale(arrival_day)         0.127      0.329    0.39    0.699  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
simulationOutput <- simulateResiduals(fittedModel = success_m3, plot=F)
plot(simulationOutput)

Similar effects when males near predated nests excluded.

success_m4_flex <- as_flextable(success_m4)
if(save_tables==TRUE){
save_as_docx(success_m4_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/odba_reproduction_model.docx"))
}
m4_effects = effects::allEffects(success_m4, xlevels = 100)
m4_active_eff = m4_effects$`scale(mean_ODBA)` %>% data.frame
m4_tenure_eff = m4_effects$`scale(log(tenure_days))` %>% data.frame

Fig: ODBA and reproductive success

egg_odba_plot <- ggplot(ind_fertile_activity)+
  geom_line(data = m4_active_eff,
            aes(y = fit, x = mean_ODBA),
            colour=col_scale[2]) +
  geom_ribbon(data = m4_active_eff,
              aes(y = fit, x = mean_ODBA, ymin = lower, ymax = upper),
              alpha = 0.2, fill=col_scale[2])+
  geom_jitter(aes(x=mean_ODBA, y=n_eggs, shape=likely_extra_nests),
             fill=male_colour,
             colour=male_colour,
             size=point_size,
             alpha=point_alpha,
             width=0.005, height=0.15) +
  theme_classic()+
  scale_shape_manual(values=c(default_shape, default_shape), guide="none")+
  scale_x_continuous(expand = c(0, 0), breaks=c(seq(0,0.5,by=0.1))) +
  coord_cartesian(xlim=c(0.1,0.52), ylim=c(0,6.2))+
  xlab("\n mean ODBA (g)") +
  ylab("reproductive success\n")+
  theme(legend.position="right",
        plot.margin=unit(c(0,0,0,0), 'lines'),
        text=element_text(family=font.family),
        legend.title = element_text(size=legend.heading.size),
        legend.text = element_text(size=label.size),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size))

egg_odba_plot

Combined figure: Mating success and behaviour

fig3 <- success_active_plot +
  success_odba_plot + ylab("\n") +
  success_tenure_plot + ylab("\n") +
  plot_layout(nrow = 1,
              widths=c(1,1.05,1.05),
              heights=c(1))+
  plot_annotation(tag_levels = 'a',
                  tag_prefix = '(',
                  tag_suffix = ')') &
  theme(plot.tag = element_text(size=label.size, family = font.family, face="italic"),
        plot.tag.position = c(.3, .96))

fig3

if(save_plots==TRUE){
ggsave(plot=fig3,
       filename='figure3-beh-mating-success.jpg',
       path=output_path,
       width=16,
       height=7,
       units="cm",
       dpi=300)
}

Combined supplementary figure: Reproductive success and behaviour

supfigsuccess <- egg_active_plot +
    egg_odba_plot + ylab("\n") +
  egg_tenure_plot + ylab("\n") +
  plot_layout(nrow = 1,
              widths=c(1,1.05,1.05),
              heights=c(1))+
  plot_annotation(tag_levels = 'a',
                  tag_prefix = '(',
                  tag_suffix = ')') &
  theme(plot.tag = element_text(size=label.size, family = font.family, face="italic"),
        plot.tag.position = c(.3, .96))

supfigsuccess

if(save_plots==TRUE){
ggsave(plot=supfigsuccess,
       filename='figure-rsuccess.jpg',
       path=output_path,
       width=16,
       height=7,
       units="cm",
       dpi=300)
}
rm(m1_active_eff, m1_active_eff2, m1_effects, m1_effects2, m1_tenure_eff,
   m2_active_eff, m2_effects, m2_tenure_eff,
   m3_active_eff, m3_effects, m3_tenure_eff,
   m4_active_eff, m4_effects, m4_tenure_eff,
   success_m1, success_m1_flex,
   success_m2, success_m2_flex,
   success_m3, success_m3_flex,
   success_m4, success_m4_flex,
   success_mt, success_mt_flex)

Q: Does activity relative to rivals predict success at a given nest?

For this analysis, we don’t want to limit the dataset just to males with at least 2 days of data.

But we do want to exclude nests that were depredated before genotyping.

We adjust the dataset accordingly.

Median daily distance from the nest site tended to be within 200 metres (in 14/18 cases).

The table below provides a summary of how many days each male who sired offspring was recorded, during the days when the female he mated with was fertile.

nest_male_summary <- nest_male_date %>%
  mutate(dist_m=as.numeric(dist_m),
         success=as.integer(success)) %>%
  filter(success==1) %>%
  select(nest, date, bird_ID, success, colour_combo) %>%
  group_by(nest, bird_ID, success, colour_combo) %>%
  summarise(n_days_father_recorded=n()) %>%
  arrange(n_days_father_recorded) %>%
  as.data.frame()

nest_male_summary <- nest_dates %>%
  group_by(nest) %>%
  summarise(n_days_mother_fertile=n()) %>%
  full_join(., nest_parentage, by='nest') %>%
  filter(father1!="unbanded") %>%
  full_join(., nest_male_summary, by='nest') %>%
  select(nest, colour_combo, n_days_mother_fertile, n_days_father_recorded) %>%
  mutate(father_unrecorded_days=n_days_mother_fertile-n_days_father_recorded) %>%
  arrange(father_unrecorded_days)

nest_male_summary %>% flextable()

nest

colour_combo

n_days_mother_fertile

n_days_father_recorded

father_unrecorded_days

P1401

Y,Y,Y

8

8

0

P203

DB,R,O

7

7

0

P214

O,DG,DG

8

8

0

P1003

O,DB,O

8

7

1

P208

O,Y,O

8

7

1

P401

PI,PI,R

8

7

1

P403

PI,R,PI

8

7

1

P1905

DG,O,O

8

6

2

P201

Y,DG,PI

8

5

3

P1903

DB,DG,Y

8

4

4

P701

PI,PI,R

8

4

4

P204

DB,R,R

8

3

5

P403

O,DG,DG

8

3

5

P702

Y,Y,R

8

2

6

P1001

DG,DB,DG

8

1

7

P1103

DG,R,R

8

1

7

P1403

Y,DG,PI

8

1

7

P701

R,DG,PI

8

1

7

nest_male_summary %>%
  summarise(mean_unrecorded_days=mean(father_unrecorded_days),
            sd_unrecorded_days=sd(father_unrecorded_days))
## # A tibble: 1 × 2
##   mean_unrecorded_days sd_unrecorded_days
##                  <dbl>              <dbl>
## 1                 3.39               2.70

Model: Activity relative to competitors

nest_male_date_sub <- nest_male_date %>%
  filter(!nest %in% predated_nests$nest) %>%
  mutate(dist_m=as.numeric(dist_m),
         success=as.integer(success)) %>%
  select(nest, date, bird_ID, colour_combo, mean_ODBA, prop_activity, dist_m, success) %>%
  filter(dist_m <= 200) %>% 
  group_by(date, nest) %>% 
  mutate(z_prop_act = scale(prop_activity),
         z_ODBA = scale(mean_ODBA)) %>%
  filter(!is.na(z_prop_act)) %>%
  as.data.frame()
nest_date_n <- nest_male_date_sub %>%
  group_by(nest, date) %>%
  summarise(n=n()) %>%
  arrange(n)

hist(nest_date_n$n, breaks=20)

Median sample size (number of rival males per nest per day):

median(nest_date_n$n)
## [1] 4

Mean sample size:

mean(nest_date_n$n)
## [1] 4.48

Standard deviation:

sd(nest_date_n$n)
## [1] 2.14
plot(nest_male_date_sub$dist_m, nest_male_date_sub$z_prop_act)

nest_male_m <- glmmTMB(success ~ scale(z_prop_act) + 
                         (1|bird_ID) +
                         (1|nest),
                       data = nest_male_date_sub,
                       family = binomial())

summary(nest_male_m)
##  Family: binomial  ( logit )
## Formula:          success ~ scale(z_prop_act) + (1 | bird_ID) + (1 | nest)
## Data: nest_male_date_sub
## 
##      AIC      BIC   logLik deviance df.resid 
##      254      273     -123      246      987 
## 
## Random effects:
## 
## Conditional model:
##  Groups  Name        Variance Std.Dev.
##  bird_ID (Intercept) 117      10.8    
##  nest    (Intercept) 114      10.7    
## Number of obs: 991, groups:  bird_ID, 62; nest, 40
## 
## Conditional model:
##                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)        -20.179      4.671   -4.32 0.000016 ***
## scale(z_prop_act)    0.734      0.311    2.36    0.018 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
simulationOutput <- simulateResiduals(fittedModel = nest_male_m, plot=F)
plot(simulationOutput)

testZeroInflation(simulationOutput)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1, p-value = 0.3
## alternative hypothesis: two.sided
testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.6, p-value = 0.2
## alternative hypothesis: two.sided
nm_rel_effects <- confint(nest_male_m)
nm_rel_effects <- cbind(name=rownames(nm_rel_effects), 
                        data.frame(nm_rel_effects,
                                   row.names=NULL)) %>%
  rename(`2.5`=`X2.5..`,
         `97.5`=`X97.5..`,
         estimate=Estimate) %>%
  mutate(effect_type="rel")
if(save_tables==TRUE){
nest_m_flex <- as_flextable(nest_male_m)
save_as_docx(nest_m_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/scaled_nest_male_model.docx"))
}

Fig: Activity relative to competitors

nest_male_plotdat <- rbind(nest_male_date_sub %>% mutate(plot_type="box"), 
                              nest_male_date_sub %>% mutate(plot_type="point"))

nest_male_plotdat$plot_factor <- factor(
  paste0("s", nest_male_plotdat$success, "-", nest_male_plotdat$plot_type),
  levels = c("s0-point", "s0-box", "s1-box", "s1-point"))

nest_male_boxplot <- ggplot() +
  theme_classic()+
    geom_point(data = nest_male_plotdat %>% 
                 filter(plot_factor == "s0-point"),
               aes(x = plot_factor,
                   y = z_prop_act),
             size=small.point.size,
             alpha=point_alpha,
             shape=default_shape,
             fill=male_colour,
             position=position_jitter(w=0.08, h=0, seed = 1),
             ) +
  geom_boxplot(data = nest_male_plotdat %>% 
                 filter(plot_factor == "s0-box"),
               aes(x = plot_factor,
                   y = z_prop_act),
               outlier.shape=default_shape,
               fill="white",
               width=0.7) +
  geom_boxplot(data = nest_male_plotdat %>% 
                 filter(plot_factor == "s1-box"),
               aes(x = plot_factor,
                   y = z_prop_act),
               outlier.shape=default_shape,
               fill="white",
               width=0.7)+
  geom_point(data = nest_male_plotdat %>% 
                 filter(plot_factor == "s1-point"),
               aes(x = plot_factor,
                   y = z_prop_act),
             size=small.point.size,
             alpha=point_alpha,
             shape=default_shape,
             fill=male_colour,
             position=position_jitter(w=0.08, h=0, seed = 1),
             ) +
  coord_cartesian(ylim=c(-2.8,2.8))+
  scale_y_continuous(breaks=c(seq(-2,2,by=1)),
                     expand=c(0,0))+
  scale_x_discrete(breaks = c("s0-point", "s1-point"),
                   labels = c("unsuccessful", "successful"))+
  xlab("\n mating success")+
  ylab("male activity\n (scaled relative to local males)\n")+
  theme(text=element_text(family=font.family),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size),
        legend.title = element_text(size=legend.heading.size),
        legend.text = element_text(size=label.size),
        legend.position="right",
        legend.justification="top",
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "grey"))

nest_male_boxplot

Model: ODBA relative to competitors

nest_male_m2 <- glmmTMB(success ~ scale(z_ODBA) + 
                         (1|bird_ID) +
                         (1|nest),
                       data = nest_male_date_sub,
                       family = binomial())

summary(nest_male_m2)
##  Family: binomial  ( logit )
## Formula:          success ~ scale(z_ODBA) + (1 | bird_ID) + (1 | nest)
## Data: nest_male_date_sub
## 
##      AIC      BIC   logLik deviance df.resid 
##      252      272     -122      244      987 
## 
## Random effects:
## 
## Conditional model:
##  Groups  Name        Variance Std.Dev.
##  bird_ID (Intercept) 120      11      
##  nest    (Intercept) 121      11      
## Number of obs: 991, groups:  bird_ID, 62; nest, 40
## 
## Conditional model:
##               Estimate Std. Error z value  Pr(>|z|)    
## (Intercept)    -20.575      4.577   -4.50 0.0000069 ***
## scale(z_ODBA)    0.912      0.353    2.58    0.0098 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
simulationOutput <- simulateResiduals(fittedModel = nest_male_m2, plot=F)
plot(simulationOutput)

testZeroInflation(simulationOutput)

## 
##  DHARMa zero-inflation test via comparison to expected zeros with
##  simulation under H0 = fitted model
## 
## data:  simulationOutput
## ratioObsSim = 1, p-value = 0.3
## alternative hypothesis: two.sided
testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 0.6, p-value = 0.3
## alternative hypothesis: two.sided
nm_rel_effects2 <- confint(nest_male_m2)
nm_rel_effects2 <- cbind(name=rownames(nm_rel_effects2), 
                        data.frame(nm_rel_effects2,
                                   row.names=NULL)) %>%
  rename(`2.5`=`X2.5..`,
         `97.5`=`X97.5..`,
         estimate=Estimate) %>%
  mutate(effect_type="rel")
if(save_tables==TRUE){
nest_m_flex <- as_flextable(nest_male_m2)
save_as_docx(nest_m_flex, 
             pr_section = sect_properties,
             path = glue("{output_path}/scaled_nest_male_odba_model.docx"))
}

Fig: ODBA relative to competitors

nest_male_boxplot2 <- ggplot() +
  theme_classic()+
    geom_point(data = nest_male_plotdat %>% 
                 filter(plot_factor == "s0-point"),
               aes(x = plot_factor,
                   y = z_ODBA),
             size=small.point.size,
             alpha=point_alpha,
             shape=default_shape,
             fill=male_colour,
             position=position_jitter(w=0.08, h=0, seed = 1),
             ) +
  geom_boxplot(data = nest_male_plotdat %>% 
                 filter(plot_factor == "s0-box"),
               aes(x = plot_factor,
                   y = z_ODBA),
               outlier.shape=default_shape,
               fill="white",
               width=0.7) +
  geom_boxplot(data = nest_male_plotdat %>% 
                 filter(plot_factor == "s1-box"),
               aes(x = plot_factor,
                   y = z_ODBA),
               outlier.shape=default_shape,
               fill="white",
               width=0.7)+
  geom_point(data = nest_male_plotdat %>% 
                 filter(plot_factor == "s1-point"),
               aes(x = plot_factor,
                   y = z_ODBA),
             size=small.point.size,
             alpha=point_alpha,
             shape=default_shape,
             fill=male_colour,
             position=position_jitter(w=0.08, h=0, seed = 1),
             ) +
  coord_cartesian(ylim=c(-2.8,2.8))+
  scale_y_continuous(breaks=c(seq(-2,2,by=1)),
                     expand=c(0,0))+
  scale_x_discrete(breaks = c("s0-point", "s1-point"),
                   labels = c("unsuccessful", "successful"))+
  xlab("\n mating success")+
  ylab("male ODBA\n (scaled relative to local males)\n")+
  theme(text=element_text(family=font.family),
        axis.title = element_text(size=axis.heading.size),
        axis.text = element_text(size=label.size),
        legend.title = element_text(size=legend.heading.size),
        legend.text = element_text(size=label.size),
        legend.position="right",
        legend.justification="top",
        legend.background = element_blank(),
        legend.box.background = element_rect(colour = "grey"))

nest_male_boxplot2

Combined figure: Activity and ODBA relative to competitors

fig4 <- nest_male_boxplot +
  plot_spacer() +
  nest_male_boxplot2 +
  plot_layout(nrow = 1,
              widths=c(1,0.1,1),
              heights=c(1))+
  plot_annotation(tag_levels = 'a',
                  tag_prefix = '(',
                  tag_suffix = ')') &
  theme(plot.tag = element_text(size=label.size, family = font.family, face="italic"),
        plot.tag.position = c(.28, .96))

fig4

if(save_plots==TRUE){
ggsave(plot=fig4,
       filename='figure4-relative-male-activity.jpg',
       path=output_path,
       width=16,
       height=7,
       units="cm",
       dpi=300)
}

Q: How repeatable is activity during the fertile period?

repeat_lme <- lmer(prop_activity ~ (1|bird_ID),
                      data = daily_activity %>% filter(breeding_stage=="fertile"))

summary(repeat_lme)
## Linear mixed model fit by REML ['lmerMod']
## Formula: prop_activity ~ (1 | bird_ID)
##    Data: daily_activity %>% filter(breeding_stage == "fertile")
## 
## REML criterion at convergence: -313
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7510 -0.5959 -0.0531  0.5941  2.5074 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  bird_ID  (Intercept) 0.0169   0.130   
##  Residual             0.0260   0.161   
## Number of obs: 555, groups:  bird_ID, 81
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   0.6055     0.0162    37.4
simulationOutput <- simulateResiduals(fittedModel = repeat_lme, plot=F)
testDispersion(simulationOutput)

## 
##  DHARMa nonparametric dispersion test via sd of residuals fitted vs.
##  simulated
## 
## data:  simulationOutput
## dispersion = 1, p-value = 0.9
## alternative hypothesis: two.sided
# Get variance
vc = VarCorr(repeat_lme)
withingroup_var=attr(vc, 'sc')[1]^2
betweengroup_var=attr(vc$bird_ID, 'stddev')[1]^2

# repeatability
R = betweengroup_var/(betweengroup_var+withingroup_var)
R
## (Intercept) 
##       0.394

Within-individual repeatability of daily activity during the fertile phase was R = 0.394.